4ad7f602af8b04c504f57f4939181edcb8c40057
[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 "AliVEvent.h"
41 #include "AliESD.h"
42 #include "AliESDEvent.h"
43 #include "AliESDHeader.h"
44 #include "AliESDInputHandler.h"
45 #include "AliESDZDC.h"
46 #include "AliESDFMD.h"
47 #include "AliESDVZERO.h"
48 #include "AliCentrality.h"
49 #include "AliESDtrackCuts.h"
50 #include "AliMultiplicity.h"
51 #include "AliAODHandler.h"
52 #include "AliAODEvent.h"
53 #include "AliESDVertex.h"
54 #include "AliAODVertex.h"
55 #include "AliAODMCHeader.h"
56 #include "AliMCEvent.h"
57 #include "AliMCEventHandler.h"
58 #include "AliMCParticle.h"
59 #include "AliStack.h"
60 #include "AliHeader.h"
61 #include "AliAODMCParticle.h"
62 #include "AliAnalysisTaskSE.h"
63 #include "AliGenEventHeader.h"
64 #include "AliGenHijingEventHeader.h"
65 #include "AliPhysicsSelectionTask.h"
66 #include "AliPhysicsSelection.h"
67 #include "AliBackgroundSelection.h"
68 #include "AliESDUtils.h"
69
70 ClassImp(AliCentralitySelectionTask)
71
72
73 //________________________________________________________________________
74 AliCentralitySelectionTask::AliCentralitySelectionTask():
75 AliAnalysisTaskSE(),
76   fDebug(0),
77   fAnalysisInput("ESD"),
78   fIsMCInput(kFALSE),
79   fFile(0),
80   fFile2(0),
81   fCurrentRun(-1),
82   fRunNo(-1),
83   fLowRunN(0),
84   fHighRunN(0),
85   fUseScaling(0),
86   fUseCleaning(0),
87   fTrackCuts(0),
88   fZVCut(10),
89   fOutliersCut(5),
90   fQuality(0),
91   fCentV0M(0),
92   fCentFMD(0),
93   fCentTRK(0),
94   fCentTKL(0),
95   fCentCL0(0),
96   fCentCL1(0),
97   fCentV0MvsFMD(0),
98   fCentTKLvsV0M(0),
99   fCentZEMvsZDC(0),
100   fHtempV0M(0),
101   fHtempFMD(0),
102   fHtempTRK(0),
103   fHtempTKL(0),
104   fHtempCL0(0),
105   fHtempCL1(0),
106   fHtempV0MvsFMD(0),
107   fHtempTKLvsV0M(0),
108   fHtempZEMvsZDC(0),
109   fOutputList(0),
110   fHOutCentV0M     (0),
111   fHOutCentFMD     (0),
112   fHOutCentTRK     (0),
113   fHOutCentTKL     (0),
114   fHOutCentCL0     (0),
115   fHOutCentCL1     (0),
116   fHOutCentV0MvsFMD(0),
117   fHOutCentTKLvsV0M(0),
118   fHOutCentZEMvsZDC(0),
119   fHOutCentV0MvsCentCL1(0),
120   fHOutCentV0MvsCentTRK(0),
121   fHOutCentTRKvsCentCL1(0),
122   fHOutMultV0M(0),
123   fHOutMultV0R(0),
124   fHOutMultFMD(0),
125   fHOutMultTRK(0),
126   fHOutMultTKL(0),
127   fHOutMultCL0(0),
128   fHOutMultCL1(0),
129   fHOutMultV0MvsZDN(0),
130   fHOutMultZEMvsZDN(0),
131   fHOutMultV0MvsZDC(0),
132   fHOutMultZEMvsZDC(0),
133   fHOutMultV0MvsCL1(0),
134   fHOutMultV0MvsTRK(0),
135   fHOutMultTRKvsCL1(0),
136   fHOutCentV0M_qual1(0),
137   fHOutCentTRK_qual1(0),
138   fHOutCentCL1_qual1(0),
139   fHOutCentV0M_qual2(0),
140   fHOutCentTRK_qual2(0),
141   fHOutCentCL1_qual2(0),
142   fHOutQuality(0),
143   fHOutVertex(0)
144 {   
145   // Default constructor
146   AliInfo("Centrality Selection enabled.");
147   fLowRunN =136851;
148   fHighRunN=139517;
149
150   for (Int_t i=0; i < 2667; i++) {
151     V0MScaleFactor[i]=0.0;
152     SPDScaleFactor[i]=0.0;
153     TPCScaleFactor[i]=0.0;
154     V0MScaleFactorMC[i]=0.0;
155   }
156   fUseScaling=kTRUE;
157   fUseCleaning=kTRUE;
158 }   
159
160 //________________________________________________________________________
161 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
162   AliAnalysisTaskSE(name),
163   fDebug(0),
164   fAnalysisInput("ESD"),
165   fIsMCInput(kFALSE),
166   fFile(0),
167   fFile2(0),
168   fCurrentRun(-1),
169   fRunNo(-1),
170   fLowRunN(0),
171   fHighRunN(0),
172   fUseScaling(0),
173   fUseCleaning(0),
174   fTrackCuts(0),
175   fZVCut(10),
176   fOutliersCut(5),
177   fQuality(0),
178   fCentV0M(0),
179   fCentFMD(0),
180   fCentTRK(0),
181   fCentTKL(0),
182   fCentCL0(0),
183   fCentCL1(0),
184   fCentV0MvsFMD(0),
185   fCentTKLvsV0M(0),
186   fCentZEMvsZDC(0),
187   fHtempV0M(0),
188   fHtempFMD(0),
189   fHtempTRK(0),
190   fHtempTKL(0),
191   fHtempCL0(0),
192   fHtempCL1(0),
193   fHtempV0MvsFMD(0),
194   fHtempTKLvsV0M(0),
195   fHtempZEMvsZDC(0),
196   fOutputList(0),
197   fHOutCentV0M     (0),
198   fHOutCentFMD     (0),
199   fHOutCentTRK     (0),
200   fHOutCentTKL     (0),
201   fHOutCentCL0     (0),
202   fHOutCentCL1     (0),
203   fHOutCentV0MvsFMD(0),
204   fHOutCentTKLvsV0M(0),
205   fHOutCentZEMvsZDC(0),
206   fHOutCentV0MvsCentCL1(0),
207   fHOutCentV0MvsCentTRK(0),
208   fHOutCentTRKvsCentCL1(0),
209   fHOutMultV0M(0),
210   fHOutMultV0R(0),
211   fHOutMultFMD(0),
212   fHOutMultTRK(0),
213   fHOutMultTKL(0),
214   fHOutMultCL0(0),
215   fHOutMultCL1(0),
216   fHOutMultV0MvsZDN(0),
217   fHOutMultZEMvsZDN(0),
218   fHOutMultV0MvsZDC(0),
219   fHOutMultZEMvsZDC(0),
220   fHOutMultV0MvsCL1(0),
221   fHOutMultV0MvsTRK(0),
222   fHOutMultTRKvsCL1(0),
223   fHOutCentV0M_qual1(0),
224   fHOutCentTRK_qual1(0),
225   fHOutCentCL1_qual1(0),
226   fHOutCentV0M_qual2(0),
227   fHOutCentTRK_qual2(0),
228   fHOutCentCL1_qual2(0),
229   fHOutQuality(0),
230   fHOutVertex(0)
231 {
232     // Default constructor
233     fLowRunN =136851;
234     fHighRunN=139517;
235
236   AliInfo("Centrality Selection enabled.");
237   DefineOutput(1, TList::Class());
238   for (Int_t i=0; i<2667; i++) {
239     V0MScaleFactor[i]=0.0;
240     SPDScaleFactor[i]=0.0;
241     TPCScaleFactor[i]=0.0;
242     V0MScaleFactorMC[i]=0.0;
243   }
244   fUseScaling=kTRUE;
245   fUseCleaning=kTRUE;
246 }
247
248 //________________________________________________________________________
249 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
250 {
251   // Assignment operator
252   if (this!=&c) {
253     AliAnalysisTaskSE::operator=(c);
254   }
255   return *this;
256 }
257
258 //________________________________________________________________________
259 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
260   AliAnalysisTaskSE(ana),
261   fDebug(ana.fDebug),     
262   fAnalysisInput(ana.fDebug),
263   fIsMCInput(ana.fIsMCInput),
264   fFile(ana.fFile),
265   fFile2(ana.fFile2),
266   fCurrentRun(ana.fCurrentRun),
267   fRunNo(ana.fRunNo),
268   fLowRunN(ana.fLowRunN),
269   fHighRunN(ana.fHighRunN),
270   fUseScaling(ana.fUseScaling),
271   fUseCleaning(ana.fUseCleaning),
272   fTrackCuts(ana.fTrackCuts),
273   fZVCut(ana.fZVCut),
274   fOutliersCut(ana.fOutliersCut),
275   fQuality(ana.fQuality),
276   fCentV0M(ana.fCentV0M),
277   fCentFMD(ana.fCentFMD),
278   fCentTRK(ana.fCentTRK),
279   fCentTKL(ana.fCentTKL),
280   fCentCL0(ana.fCentCL0),
281   fCentCL1(ana.fCentCL1),
282   fCentV0MvsFMD(ana.fCentV0MvsFMD),
283   fCentTKLvsV0M(ana.fCentTKLvsV0M),
284   fCentZEMvsZDC(ana.fCentZEMvsZDC),
285   fHtempV0M(ana.fHtempV0M),
286   fHtempFMD(ana.fHtempFMD),
287   fHtempTRK(ana.fHtempTRK),
288   fHtempTKL(ana.fHtempTKL),
289   fHtempCL0(ana.fHtempCL0),
290   fHtempCL1(ana.fHtempCL1),
291   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
292   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
293   fHtempZEMvsZDC(ana.fHtempZEMvsZDC),
294   fOutputList(ana.fOutputList),
295   fHOutCentV0M     (ana.fHOutCentV0M     ),
296   fHOutCentFMD     (ana.fHOutCentFMD     ),
297   fHOutCentTRK     (ana.fHOutCentTRK     ),
298   fHOutCentTKL     (ana.fHOutCentTKL     ),
299   fHOutCentCL0     (ana.fHOutCentCL0     ),
300   fHOutCentCL1     (ana.fHOutCentCL1     ),
301   fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD),
302   fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M),
303   fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC),
304   fHOutCentV0MvsCentCL1(ana.fHOutCentV0MvsCentCL1),
305   fHOutCentV0MvsCentTRK(ana.fHOutCentV0MvsCentTRK),
306   fHOutCentTRKvsCentCL1(ana.fHOutCentTRKvsCentCL1),
307   fHOutMultV0M(ana.fHOutMultV0M),
308   fHOutMultV0R(ana.fHOutMultV0R),
309   fHOutMultFMD(ana.fHOutMultFMD),
310   fHOutMultTRK(ana.fHOutMultTRK),
311   fHOutMultTKL(ana.fHOutMultTKL),
312   fHOutMultCL0(ana.fHOutMultCL0),
313   fHOutMultCL1(ana.fHOutMultCL1),
314   fHOutMultV0MvsZDN(ana.fHOutMultV0MvsZDN),
315   fHOutMultZEMvsZDN(ana.fHOutMultZEMvsZDN),
316   fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC),
317   fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),
318   fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),
319   fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),
320   fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1),
321   fHOutCentV0M_qual1(ana.fHOutCentV0M_qual1),
322   fHOutCentTRK_qual1(ana.fHOutCentTRK_qual1),
323   fHOutCentCL1_qual1(ana.fHOutCentCL1_qual1),
324   fHOutCentV0M_qual2(ana.fHOutCentV0M_qual2),
325   fHOutCentTRK_qual2(ana.fHOutCentTRK_qual2),
326   fHOutCentCL1_qual2(ana.fHOutCentCL1_qual2),
327   fHOutQuality(ana.fHOutQuality),
328   fHOutVertex(ana.fHOutVertex)
329 {
330   // Copy Constructor   
331     for (Int_t i=0; i<2667; i++) {
332         V0MScaleFactor[i]=0.0;
333         SPDScaleFactor[i]=0.0;
334         TPCScaleFactor[i]=0.0;
335         V0MScaleFactorMC[i]=0.0;
336     }
337
338 }
339  
340 //________________________________________________________________________
341 AliCentralitySelectionTask::~AliCentralitySelectionTask()
342 {
343   // Destructor  
344   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
345   if (fTrackCuts) delete fTrackCuts;
346 }  
347
348 //________________________________________________________________________
349 void AliCentralitySelectionTask::UserCreateOutputObjects()
350 {  
351   // Create the output containers
352   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
353   AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
354
355   fOutputList = new TList();
356   fOutputList->SetOwner();
357   fHOutCentV0M     = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",501,0,101);
358   fHOutCentFMD     = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",501,0,101);
359   fHOutCentTRK     = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",501,0,101);
360   fHOutCentTKL     = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",501,0,101);
361   fHOutCentCL0     = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",501,0,101);
362   fHOutCentCL1     = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",501,0,101);
363   fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",501,0,101);
364   fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",501,0,101);
365   fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",501,0,101);
366   fHOutCentV0MvsCentCL1= new TH2F("fHOutCentV0MvsCentCL1","fHOutCentV0MvsCentCL1; Cent V0 vs Cent SPD",501,0,101,501,0,101);
367   fHOutCentV0MvsCentTRK= new TH2F("fHOutCentV0MvsCentTRK","fHOutCentV0MvsCentTRK; Cent V0 vs Cent TPC",501,0,101,501,0,101);
368   fHOutCentTRKvsCentCL1= new TH2F("fHOutCentTRKvsCentCL1","fHOutCentTRKvsCentCL1; Cent TPC vs Cent SPD",501,0,101,501,0,101);
369
370   fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,25000);
371   fHOutMultV0R = new TH1F("fHOutMultV0R","fHOutMultV0R; Multiplicity V0",30000,0,30000);
372   fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000);
373   fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000);
374   fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000);
375   fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000);
376   fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000);
377   fHOutMultV0MvsZDN = new TH2F("fHOutMultV0MvsZDN","fHOutMultV0MvsZDN; Multiplicity V0; Energy ZDC-N",500,0,25000,500,0,180000);
378   fHOutMultZEMvsZDN = new TH2F("fHOutMultZEMvsZDN","fHOutMultZEMvsZDN; Energy ZEM; Energy ZDC-N",500,0,2500,500,0,180000);
379   fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,25000,500,0,200000);
380   fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,200000);
381   fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
382   fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
383   fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
384
385   fHOutCentV0M_qual1 = new TH1F("fHOutCentV0M_qual1","fHOutCentV0M_qual1; Centrality V0",501,0,100);
386   fHOutCentTRK_qual1 = new TH1F("fHOutCentTRK_qual1","fHOutCentTRK_qual1; Centrality TPC",501,0,100);
387   fHOutCentCL1_qual1 = new TH1F("fHOutCentCL1_qual1","fHOutCentCL1_qual1; Centrality SPD outer",501,0,100);
388
389   fHOutCentV0M_qual2 = new TH1F("fHOutCentV0M_qual2","fHOutCentV0M_qual2; Centrality V0",101,-0.1,100.1);
390   fHOutCentTRK_qual2 = new TH1F("fHOutCentTRK_qual2","fHOutCentTRK_qual2; Centrality TPC",101,-0.1,100.1);
391   fHOutCentCL1_qual2 = new TH1F("fHOutCentCL1_qual2","fHOutCentCL1_qual2; Centrality SPD outer",101,-0.1,100.1);
392
393   fHOutQuality = new TH1F("fHOutQuality", "fHOutQuality", 10,-0.5,9.5);
394   fHOutVertex  = new TH1F("fHOutVertex", "fHOutVertex", 100,-20,20);
395
396   fOutputList->Add(  fHOutCentV0M     );
397   fOutputList->Add(  fHOutCentFMD     );
398   fOutputList->Add(  fHOutCentTRK     );
399   fOutputList->Add(  fHOutCentTKL     );
400   fOutputList->Add(  fHOutCentCL0     );
401   fOutputList->Add(  fHOutCentCL1     );
402   fOutputList->Add(  fHOutCentV0MvsFMD);
403   fOutputList->Add(  fHOutCentTKLvsV0M);
404   fOutputList->Add(  fHOutCentZEMvsZDC);
405   fOutputList->Add(  fHOutCentV0MvsCentCL1);
406   fOutputList->Add(  fHOutCentV0MvsCentTRK);
407   fOutputList->Add(  fHOutCentTRKvsCentCL1);
408   fOutputList->Add(  fHOutMultV0M); 
409   fOutputList->Add(  fHOutMultV0R); 
410   fOutputList->Add(  fHOutMultFMD); 
411   fOutputList->Add(  fHOutMultTRK); 
412   fOutputList->Add(  fHOutMultTKL); 
413   fOutputList->Add(  fHOutMultCL0); 
414   fOutputList->Add(  fHOutMultCL1); 
415   fOutputList->Add(  fHOutMultV0MvsZDN);
416   fOutputList->Add(  fHOutMultZEMvsZDN);
417   fOutputList->Add(  fHOutMultV0MvsZDC);
418   fOutputList->Add(  fHOutMultZEMvsZDC);
419   fOutputList->Add(  fHOutMultV0MvsCL1);
420   fOutputList->Add(  fHOutMultV0MvsTRK);
421   fOutputList->Add(  fHOutMultTRKvsCL1);
422   fOutputList->Add(  fHOutCentV0M_qual1 );
423   fOutputList->Add(  fHOutCentTRK_qual1 );
424   fOutputList->Add(  fHOutCentCL1_qual1 );                   
425   fOutputList->Add(  fHOutCentV0M_qual2 );
426   fOutputList->Add(  fHOutCentTRK_qual2 );
427   fOutputList->Add(  fHOutCentCL1_qual2 );
428   fOutputList->Add(  fHOutQuality );
429   fOutputList->Add(  fHOutVertex );
430
431
432   fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
433
434   PostData(1, fOutputList); 
435
436   MyInitScaleFactor();
437   if (fIsMCInput) MyInitScaleFactorMC();
438 }
439
440 //________________________________________________________________________
441 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
442
443   // Execute analysis for current event:
444   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
445   
446   Float_t  zncEnergy = 0.;          //  ZNC Energy
447   Float_t  zpcEnergy = 0.;          //  ZPC Energy
448   Float_t  znaEnergy = 0.;          //  ZNA Energy
449   Float_t  zpaEnergy = 0.;          //  ZPA Energy
450   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
451   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
452   Bool_t   zdcEnergyCal = kFALSE;   // if zdc is calibrated (in pass2)
453
454   Int_t    nTracks = 0;             //  no. tracks
455   Int_t    nTracklets = 0;          //  no. tracklets
456   Int_t    nClusters[6] = {0};      //  no. clusters on 6 ITS layers
457   Int_t    nChips[2];               //  no. chips on 2 SPD layers
458   Float_t  spdCorr =0;              //  corrected spd2 multiplicity
459
460   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
461   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
462   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
463   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
464
465   Short_t v0Corr = 0;               // corrected V0 multiplicity
466   Short_t v0CorrResc = 0;           // corrected and rescaled V0 multiplicity
467
468   Float_t zvtx =0;                  // z-vertex SPD
469  
470   AliCentrality *esdCent = 0;
471
472   if(fAnalysisInput.CompareTo("ESD")==0){
473
474     AliVEvent* event = InputEvent();
475     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
476     if (!esd) {
477         AliError("No ESD Event");
478         return;
479     }
480     
481     if (fRunNo<=0) {
482       if (SetupRun(esd)<0)
483          AliFatal("Centrality File not available for this run");
484     }
485
486     esdCent = esd->GetCentrality();
487
488     // ***** V0 info    
489     AliESDVZERO* esdV0 = esd->GetVZEROData();
490     multV0A=esdV0->GetMTotV0A();
491     multV0C=esdV0->GetMTotV0C();
492
493     Float_t v0CorrR;
494     v0Corr = (Short_t)AliESDUtils::GetCorrV0(esd,v0CorrR);
495     v0CorrResc = (Short_t)v0CorrR;
496
497     // ***** Vertex Info
498     const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
499     zvtx        = vtxESD->GetZ(); 
500
501     // ***** CB info (tracklets, clusters, chips)
502     //nTracks    = event->GetNumberOfTracks();     
503     nTracks    = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
504
505     const AliMultiplicity *mult = esd->GetMultiplicity();
506
507     nTracklets = mult->GetNumberOfTracklets();
508
509     for(Int_t ilay=0; ilay<6; ilay++){
510       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
511     }
512
513     for(Int_t ilay=0; ilay<2; ilay++){
514       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
515     }
516
517     spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);    
518
519     // ***** FMD info
520     AliESDFMD *fmd = esd->GetFMDData();
521     Float_t totalMultA = 0;
522     Float_t totalMultC = 0;
523     const Float_t fFMDLowCut = 0.4;
524     
525     for(UShort_t det=1;det<=3;det++) {
526       Int_t nRings = (det==1 ? 1 : 2);
527       for (UShort_t ir = 0; ir < nRings; ir++) {          
528         Char_t   ring = (ir == 0 ? 'I' : 'O');
529         UShort_t nsec = (ir == 0 ? 20  : 40);
530         UShort_t nstr = (ir == 0 ? 512 : 256);
531         for(UShort_t sec =0; sec < nsec;  sec++)  {
532           for(UShort_t strip = 0; strip < nstr; strip++) {
533             
534             Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
535             if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
536             
537             Float_t nParticles=0;
538             
539             if(FMDmult > fFMDLowCut) {
540               nParticles = 1.;
541             }
542             
543             if (det<3) totalMultA = totalMultA + nParticles;
544             else totalMultC = totalMultC + nParticles;
545             
546           }
547         }
548       }
549     }
550     multFMDA = totalMultA;
551     multFMDC = totalMultC;
552     
553     // ***** ZDC info
554     AliESDZDC *esdZDC = esd->GetESDZDC();
555     zdcEnergyCal = esdZDC->AliESDZDC::TestBit(AliESDZDC::kEnergyCalibratedSignal);
556     if (zdcEnergyCal) {      
557       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
558       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
559       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
560       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
561     } else {
562       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
563       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
564       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
565       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
566     }
567     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
568     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
569   
570   }   
571   else if(fAnalysisInput.CompareTo("AOD")==0){
572     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
573     // to be implemented
574     printf("  AOD analysis not yet implemented!!!\n\n");
575     return;
576   } 
577
578
579   // ***** Scaling
580   // ***** Scaling for MC
581   if (fIsMCInput) {
582     fUseScaling=kFALSE;
583     Float_t temp_scalefactorV0M = MyGetScaleFactorMC(fCurrentRun);
584     v0Corr  = Short_t((multV0A+multV0C)  * temp_scalefactorV0M);
585   }
586   // ***** Scaling for Data
587   if (fUseScaling) {
588     Float_t temp_scalefactorV0M = MyGetScaleFactor(fCurrentRun,0);
589     Float_t temp_scalefactorSPD = MyGetScaleFactor(fCurrentRun,1);
590     Float_t temp_scalefactorTPC = MyGetScaleFactor(fCurrentRun,2);
591     v0Corr  = Short_t(v0Corr / temp_scalefactorV0M);
592     spdCorr = spdCorr / temp_scalefactorSPD;
593     nTracks = Int_t(nTracks / temp_scalefactorTPC);
594   }
595
596   // ***** Centrality Selection
597   if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
598   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
599   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
600   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
601   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
602   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
603   
604   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
605   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
606   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin(zem1Energy+zem2Energy,zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
607
608   // ***** Cleaning
609   if (fUseCleaning) {
610       fQuality=0;
611       fZVCut=10;
612       fOutliersCut=6;
613       
614       // ***** vertex
615       if (TMath::Abs(zvtx)>fZVCut) fQuality += 1;   
616
617       // ***** outliers
618       // **** V0 vs SPD
619       if (IsOutlierV0MSPD(spdCorr, v0Corr, int(fCentV0M))) fQuality  += 2;
620       // ***** V0 vs TPC
621       if (IsOutlierV0MTPC(nTracks, v0Corr, int(fCentV0M))) fQuality  += 4;
622       // ***** V0 vs ZDC
623        if (IsOutlierV0MZDC((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
624            (zdcEnergyCal==kFALSE) && !(fIsMCInput)) fQuality  += 8;
625        if (IsOutlierV0MZDCECal((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
626            ((zdcEnergyCal==kTRUE) || (fIsMCInput))) fQuality  += 8;
627   } else {
628       fQuality = 0;
629   }
630
631         
632   if (esdCent) {
633       esdCent->SetQuality(fQuality);
634       esdCent->SetCentralityV0M(fCentV0M);
635       esdCent->SetCentralityFMD(fCentFMD);
636       esdCent->SetCentralityTRK(fCentTRK);
637       esdCent->SetCentralityTKL(fCentTKL);
638       esdCent->SetCentralityCL0(fCentCL0);
639       esdCent->SetCentralityCL1(fCentCL1);
640       esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
641       esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
642       esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
643   }
644
645   fHOutQuality->Fill(fQuality);
646   fHOutVertex->Fill(zvtx);
647
648   fHOutMultV0MvsZDN->Fill(v0Corr,(zncEnergy+znaEnergy));
649   fHOutMultZEMvsZDN->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy));
650   fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
651   fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
652   fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
653   fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
654   
655   if (fQuality==0) {  
656     fHOutCentV0M->Fill(fCentV0M);
657     fHOutCentFMD->Fill(fCentFMD);
658     fHOutCentTRK->Fill(fCentTRK);
659     fHOutCentTKL->Fill(fCentTKL);
660     fHOutCentCL0->Fill(fCentCL0);
661     fHOutCentCL1->Fill(fCentCL1);
662     fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
663     fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
664     fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
665     fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentCL1);
666     fHOutCentV0MvsCentTRK->Fill(fCentV0M,fCentTRK);
667     fHOutCentTRKvsCentCL1->Fill(fCentTRK,fCentCL1);
668     fHOutMultV0M->Fill(v0Corr);
669     fHOutMultV0R->Fill(multV0A+multV0C);
670     fHOutMultFMD->Fill((multFMDA+multFMDC));
671     fHOutMultTRK->Fill(nTracks);
672     fHOutMultTKL->Fill(nTracklets);
673     fHOutMultCL0->Fill(nClusters[0]);
674     fHOutMultCL1->Fill(spdCorr);
675     fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
676   } else if (fQuality ==1) {
677     fHOutCentV0M_qual1->Fill(fCentV0M);
678     fHOutCentTRK_qual1->Fill(fCentTRK);
679     fHOutCentCL1_qual1->Fill(fCentCL1);
680   } else {
681     fHOutCentV0M_qual2->Fill(fCentV0M);
682     fHOutCentTRK_qual2->Fill(fCentTRK);
683     fHOutCentCL1_qual2->Fill(fCentCL1);
684   }
685
686   PostData(1, fOutputList); 
687 }
688
689 //________________________________________________________________________
690 void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename) 
691 {
692   //  Read centrality histograms
693   TDirectory *owd = gDirectory;
694   // Check if the file is present
695   TString path = gSystem->ExpandPathName(fCentfilename.Data());
696   if (gSystem->AccessPathName(path)) {
697      AliError(Form("File %s does not exist", path.Data()));
698      return;
699   }
700   fFile  = TFile::Open(fCentfilename);
701   owd->cd();
702   fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));
703   fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));
704   fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));
705   fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));
706   fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));
707   fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));
708
709   if (!fHtempV0M) AliWarning(Form("Calibration for V0M does not exist in %s", path.Data()));
710   if (!fHtempFMD) AliWarning(Form("Calibration for FMD does not exist in %s", path.Data()));
711   if (!fHtempTRK) AliWarning(Form("Calibration for TRK does not exist in %s", path.Data()));
712   if (!fHtempTKL) AliWarning(Form("Calibration for TKL does not exist in %s", path.Data()));
713   if (!fHtempCL0) AliWarning(Form("Calibration for CL0 does not exist in %s", path.Data()));
714   if (!fHtempCL1) AliWarning(Form("Calibration for CL1 does not exist in %s", path.Data()));
715   
716   owd->cd();
717 }  
718
719 //________________________________________________________________________
720 void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2) 
721 {
722   //  Read centrality histograms
723   TDirectory *owd = gDirectory;
724   TString path = gSystem->ExpandPathName(fCentfilename2.Data());
725   if (gSystem->AccessPathName(path)) {
726      AliError(Form("File %s does not exist", path.Data()));
727      return;
728   }   
729   fFile2  = TFile::Open(fCentfilename2);
730   owd->cd();
731   fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
732   fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
733   fHtempZEMvsZDC  = (TH2F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
734
735   if (!fHtempV0MvsFMD) AliWarning(Form("Calibration for V0MvsFMD does not exist in %s", path.Data()));
736   if (!fHtempTKLvsV0M) AliWarning(Form("Calibration for TKLvsV0M does not exist in %s", path.Data()));
737   if (!fHtempZEMvsZDC) AliWarning(Form("Calibration for ZEMvsZDC does not exist in %s", path.Data()));
738   
739   owd->cd();
740 }
741
742 //________________________________________________________________________
743 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
744 {
745   // Terminate analysis
746   if (fFile && fFile->IsOpen())
747     fFile->Close();  
748   if (fFile2 && fFile2->IsOpen())
749     fFile2->Close();  
750 }
751 //________________________________________________________________________
752 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* const esd)
753 {
754   // Setup files for run
755
756   if (!esd)
757     return -1;
758
759   // check if something to be done
760   if (fCurrentRun == esd->GetRunNumber())
761     return 0;
762   else
763     fCurrentRun = esd->GetRunNumber();
764   
765   AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
766
767   // CHANGE HERE FOR RUN RANGES
768   if ( fCurrentRun <= 137165 ) fRunNo = 137161;
769   else fRunNo = 137366;
770   // CHANGE HERE FOR RUN RANGES
771   
772   TString fileName(Form("%s/COMMON/CENTRALITY/data/AliCentralityBy1D_%d.root", AliAnalysisManager::GetOADBPath(), fRunNo));
773   TString fileName2(Form("%s/COMMON/CENTRALITY/data/AliCentralityByFunction_%d.root", AliAnalysisManager::GetOADBPath(), fRunNo));
774   
775   AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));
776   ReadCentralityHistos(fileName.Data());
777   ReadCentralityHistos2(fileName2.Data());
778   if (!fFile && !fFile2) {
779      AliFatal(Form("Run %d not known to centrality selection!", fCurrentRun));       
780      return -1;
781   }   
782   return 0;
783 }
784
785 //________________________________________________________________________
786 Bool_t AliCentralitySelectionTask::IsOutlierV0MSPD(Float_t spd, Float_t v0, Int_t cent) const
787 {
788 // Clean outliers
789   Float_t val= -0.143789 + 0.288874 * v0;
790   Float_t SPDsigma[101]={231.483, 189.446, 183.359, 179.923, 174.229, 170.309, 165.021, 
791                          160.84, 159.33, 154.453, 151.644, 148.337, 145.215, 142.353, 
792                          139.351, 136, 133.838, 129.885, 127.36, 125.032, 122.21, 120.3, 
793                          117.766, 114.77, 113.1, 110.268, 107.463, 105.293, 102.845, 
794                          100.835, 98.9632, 97.3287, 93.6887, 92.1066, 89.3224, 87.8382, 
795                          86.04, 83.6431, 81.9655, 80.0491, 77.8208, 76.4716, 74.2165, 
796                          72.2752, 70.4875, 68.9414, 66.8622, 65.022, 63.5134, 61.8228, 
797                          59.7166, 58.5008, 56.2789, 54.9615, 53.386, 51.2165, 49.4842, 
798                          48.259, 47.1129, 45.3115, 43.8486, 41.9207, 40.5754, 39.3872, 
799                          38.1897, 36.5401, 35.1283, 33.9702, 32.6429, 31.3612, 29.5876, 
800                          28.9319, 27.564, 26.0443, 25.2836, 23.9753, 22.8936, 21.5665, 
801                          20.7048, 19.8016, 18.7095, 18.1144, 17.2095, 16.602, 16.3233, 
802                          15.7185, 15.3006, 14.7432, 14.4174, 14.0805, 13.7638, 13.7638, 
803                          13.7638, 13.7638, 13.7638, 13.7638, 13.7638, 13.7638, 13.7638, 18.0803, 18.0803};
804
805   if ( TMath::Abs(spd-val) > fOutliersCut*SPDsigma[cent] ) 
806     return kTRUE;
807   else 
808     return kFALSE;
809 }
810
811 //________________________________________________________________________
812 Bool_t AliCentralitySelectionTask::IsOutlierV0MTPC(Int_t tracks, Float_t v0, Int_t cent) const
813 {
814 // Clean outliers
815   Float_t val = -0.540691 + 0.128358 * v0;
816   Float_t TPCsigma[101]={106.439, 89.2834, 86.7568, 85.3641, 83.379, 81.6093, 79.3189, 
817                          78.0616, 77.2167, 75.0021, 73.9957, 72.0926, 71.0442, 69.8395, 
818                          68.1169, 66.6676, 66.0038, 64.2284, 63.3845, 61.7439, 60.642, 
819                          59.5383, 58.3696, 57.0227, 56.0619, 54.7108, 53.8382, 52.3398, 
820                          51.5297, 49.9488, 49.1126, 48.208, 46.8566, 45.7724, 44.7829, 
821                          43.8726, 42.7499, 41.9307, 40.6874, 39.9619, 38.5534, 38.0884, 
822                          36.6141, 35.7482, 34.8804, 34.1769, 33.1278, 32.3435, 31.4783, 
823                          30.2587, 29.4741, 28.8575, 27.9298, 26.9752, 26.1675, 25.1234, 
824                          24.4702, 23.6843, 22.9764, 21.8579, 21.2924, 20.3241, 19.8296, 
825                          19.2465, 18.4474, 17.7216, 16.8956, 16.342, 15.626, 15.0329, 
826                          14.3911, 13.9301, 13.254, 12.6745, 12.2436, 11.7776, 11.1795, 
827                          10.673, 10.27, 9.95646, 9.50939, 9.26162, 8.95315, 8.73439, 
828                          8.67375, 8.43029, 8.34818, 8.33484, 8.40709, 8.3974, 8.32814, 
829                          8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 12.351, 12.351};
830
831   if ( TMath::Abs(tracks-val) > fOutliersCut*TPCsigma[cent] ) 
832     return kTRUE;
833   else 
834     return kFALSE;
835 }
836
837 //________________________________________________________________________
838 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDC(Float_t zdc, Float_t v0) const
839 {
840 // Clean outliers
841   Float_t val1 = 6350. - 0.26 * v0;
842   Float_t val2 = 5580.;
843   if ((zdc >  val1) || (zdc > val2)) 
844     return kTRUE;
845   else 
846     return kFALSE;
847 }
848
849 //________________________________________________________________________
850 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDCECal(Float_t /*zdc*/, Float_t /*v0*/) const
851 {
852 // Clean outliers
853     return kFALSE;
854 }
855
856 //________________________________________________________________________  
857 Float_t AliCentralitySelectionTask::MyGetScaleFactor(Int_t runnumber, Int_t flag) const
858 {
859 // Get scaling factor
860   if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
861     cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
862     return 0.0;
863   }
864
865   Float_t scalefactor=0.0;
866   if (flag==0)
867     scalefactor = V0MScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
868   else if (flag==1)
869     scalefactor = SPDScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
870   else if (flag==2)
871     scalefactor = TPCScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
872
873   return scalefactor;
874
875 }
876
877 //________________________________________________________________________  
878 Float_t AliCentralitySelectionTask::MyGetScaleFactorMC(Int_t runnumber) const
879 {
880 // Get MC scaling factor
881   if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
882     cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
883     return 0.0;
884   }
885
886   Float_t scalefactor= V0MScaleFactorMC[runnumber - fLowRunN]; // subtracting reference offset index
887   return scalefactor;
888
889 }
890
891 //________________________________________________________________________  
892 void AliCentralitySelectionTask::MyInitScaleFactor () 
893 {
894 // Initialize the scaling factors
895   for (int i=0; i<(fHighRunN-fLowRunN); i++) V0MScaleFactor[i] = 0.0;
896   for (int i=0; i<(fHighRunN-fLowRunN); i++) SPDScaleFactor[i] = 0.0;
897   for (int i=0; i<(fHighRunN-fLowRunN); i++) TPCScaleFactor[i] = 0.0;
898   
899   // scale factors determined from <V0 charge> on a run-by-run basis
900    V0MScaleFactor[310] = 1.;
901    V0MScaleFactor[311] = 1.;
902    V0MScaleFactor[514] = 1.0046;
903    V0MScaleFactor[515] = 0.983535;
904    V0MScaleFactor[579] = 0.988185;
905    V0MScaleFactor[580] = 0.983351;
906    V0MScaleFactor[581] = 0.989013;
907    V0MScaleFactor[583] = 0.990056;
908    V0MScaleFactor[588] = 0.974438;
909    V0MScaleFactor[589] = 0.981572;
910    V0MScaleFactor[590] = 0.989316;
911    V0MScaleFactor[592] = 0.98345;
912    V0MScaleFactor[688] = 0.993647;
913    V0MScaleFactor[690] = 0.994758;
914    V0MScaleFactor[693] = 0.989569;
915    V0MScaleFactor[698] = 0.993119;
916    V0MScaleFactor[744] = 0.989583;
917    V0MScaleFactor[757] = 0.990377;
918    V0MScaleFactor[787] = 0.990176;
919    V0MScaleFactor[788] = 0.98723;
920    V0MScaleFactor[834] = 1.00403;
921    V0MScaleFactor[835] = 0.994376;
922    V0MScaleFactor[840] = 0.99759;
923    V0MScaleFactor[842] = 1.01031;
924    V0MScaleFactor[900] = 0.996216;
925    V0MScaleFactor[901] = 0.994205;
926    V0MScaleFactor[992] = 0.998479;
927    V0MScaleFactor[997] = 1.00078;
928    V0MScaleFactor[1299] = 1.00515;
929    V0MScaleFactor[1303] = 1.00094;
930    V0MScaleFactor[1339] = 0.986596;
931    V0MScaleFactor[1346] = 0.972226;
932    V0MScaleFactor[1349] = 0.960358;
933    V0MScaleFactor[1350] = 0.970023;
934    V0MScaleFactor[1374] = 1.00575;
935    V0MScaleFactor[1545] = 1.00471;
936    V0MScaleFactor[1587] = 1.00611;
937    V0MScaleFactor[1588] = 1.00976;
938    V0MScaleFactor[1618] = 1.00771;
939    V0MScaleFactor[1682] = 1.01622;
940    V0MScaleFactor[1727] = 1.01305;
941    V0MScaleFactor[1728] = 1.00388;
942    V0MScaleFactor[1731] = 1.00673;
943    V0MScaleFactor[1732] = 1.00916;
944    V0MScaleFactor[1770] = 1.0092;
945    V0MScaleFactor[1773] = 1.00728;
946    V0MScaleFactor[1786] = 1.01655;
947    V0MScaleFactor[1787] = 1.00672;
948    V0MScaleFactor[1801] = 0.983339;
949    V0MScaleFactor[1802] = 1.00754;
950    V0MScaleFactor[1811] = 1.00608;
951    V0MScaleFactor[1815] = 1.01227;
952    V0MScaleFactor[1879] = 0.99907;
953    V0MScaleFactor[1881] = 0.995696;
954    V0MScaleFactor[2186] = 1.00559;
955    V0MScaleFactor[2187] = 1.00631;
956    V0MScaleFactor[2254] = 1.01512;
957    V0MScaleFactor[2256] = 0.998727;
958    V0MScaleFactor[2321] = 1.00701;
959    SPDScaleFactor[310] = 1.00211;
960    SPDScaleFactor[311] = 1.00067;
961    SPDScaleFactor[514] = 1.02268;
962    SPDScaleFactor[515] = 0.994902;
963    SPDScaleFactor[579] = 1.00215;
964    SPDScaleFactor[580] = 0.993421;
965    SPDScaleFactor[581] = 1.00129;
966    SPDScaleFactor[583] = 1.00242;
967    SPDScaleFactor[588] = 0.984762;
968    SPDScaleFactor[589] = 0.994355;
969    SPDScaleFactor[590] = 1.00073;
970    SPDScaleFactor[592] = 0.995889;
971    SPDScaleFactor[688] = 0.994532;
972    SPDScaleFactor[690] = 0.998307;
973    SPDScaleFactor[693] = 0.994052;
974    SPDScaleFactor[698] = 0.993224;
975    SPDScaleFactor[744] = 0.993279;
976    SPDScaleFactor[757] = 0.992494;
977    SPDScaleFactor[787] = 0.992678;
978    SPDScaleFactor[788] = 0.996563;
979    SPDScaleFactor[834] = 1.01116;
980    SPDScaleFactor[835] = 0.993108;
981    SPDScaleFactor[840] = 0.997574;
982    SPDScaleFactor[842] = 1.01829;
983    SPDScaleFactor[900] = 0.999438;
984    SPDScaleFactor[901] = 0.995849;
985    SPDScaleFactor[992] = 0.999227;
986    SPDScaleFactor[997] = 1.00575;
987    SPDScaleFactor[1299] = 0.99877;
988    SPDScaleFactor[1303] = 0.999682;
989    SPDScaleFactor[1339] = 0.978198;
990    SPDScaleFactor[1346] = 0.964178;
991    SPDScaleFactor[1349] = 0.959439;
992    SPDScaleFactor[1350] = 0.956945;
993    SPDScaleFactor[1374] = 0.994434;
994    SPDScaleFactor[1545] = 1.0016;
995    SPDScaleFactor[1587] = 1.00153;
996    SPDScaleFactor[1588] = 1.00698;
997    SPDScaleFactor[1618] = 1.00554;
998    SPDScaleFactor[1682] = 1.0123;
999    SPDScaleFactor[1727] = 1.011;
1000    SPDScaleFactor[1728] = 1.00143;
1001    SPDScaleFactor[1731] = 1.00486;
1002    SPDScaleFactor[1732] = 1.00646;
1003    SPDScaleFactor[1770] = 1.00515;
1004    SPDScaleFactor[1773] = 1.00485;
1005    SPDScaleFactor[1786] = 1.01215;
1006    SPDScaleFactor[1787] = 1.00419;
1007    SPDScaleFactor[1801] = 0.983327;
1008    SPDScaleFactor[1802] = 1.00529;
1009    SPDScaleFactor[1811] = 1.00367;
1010    SPDScaleFactor[1815] = 1.01045;
1011    SPDScaleFactor[1879] = 0.996374;
1012    SPDScaleFactor[1881] = 0.988827;
1013    SPDScaleFactor[2186] = 1.00354;
1014    SPDScaleFactor[2187] = 1.00397;
1015    SPDScaleFactor[2254] = 1.01138;
1016    SPDScaleFactor[2256] = 0.996641;
1017    SPDScaleFactor[2321] = 1.00357;
1018    TPCScaleFactor[310] = 1.00434;
1019    TPCScaleFactor[311] = 1.0056;
1020    TPCScaleFactor[514] = 1.02185;
1021    TPCScaleFactor[515] = 1.0036;
1022    TPCScaleFactor[579] = 1.00607;
1023    TPCScaleFactor[580] = 1.00183;
1024    TPCScaleFactor[581] = 1.00693;
1025    TPCScaleFactor[583] = 1.00746;
1026    TPCScaleFactor[588] = 0.990524;
1027    TPCScaleFactor[589] = 0.998582;
1028    TPCScaleFactor[590] = 1.00618;
1029    TPCScaleFactor[592] = 1.00088;
1030    TPCScaleFactor[688] = 1.00598;
1031    TPCScaleFactor[690] = 1.00658;
1032    TPCScaleFactor[693] = 1.00144;
1033    TPCScaleFactor[698] = 1.00279;
1034    TPCScaleFactor[744] = 1.00122;
1035    TPCScaleFactor[757] = 1.002;
1036    TPCScaleFactor[787] = 0.997818;
1037    TPCScaleFactor[788] = 0.994583;
1038    TPCScaleFactor[834] = 1.01508;
1039    TPCScaleFactor[835] = 1.00218;
1040    TPCScaleFactor[840] = 1.00569;
1041    TPCScaleFactor[842] = 1.01789;
1042    TPCScaleFactor[900] = 1.00739;
1043    TPCScaleFactor[901] = 1.00462;
1044    TPCScaleFactor[992] = 1.00502;
1045    TPCScaleFactor[997] = 1.00943;
1046    TPCScaleFactor[1299] = 1.00438;
1047    TPCScaleFactor[1303] = 0.996701;
1048    TPCScaleFactor[1339] = 0.978641;
1049    TPCScaleFactor[1346] = 0.968906;
1050    TPCScaleFactor[1349] = 0.954311;
1051    TPCScaleFactor[1350] = 0.958764;
1052    TPCScaleFactor[1374] = 0.997899;
1053    TPCScaleFactor[1545] = 0.992;
1054    TPCScaleFactor[1587] = 0.992635;
1055    TPCScaleFactor[1588] = 1.00207;
1056    TPCScaleFactor[1618] = 1.00126;
1057    TPCScaleFactor[1682] = 1.00324;
1058    TPCScaleFactor[1727] = 1.00042;
1059    TPCScaleFactor[1728] = 0.978881;
1060    TPCScaleFactor[1731] = 0.999818;
1061    TPCScaleFactor[1732] = 1.00109;
1062    TPCScaleFactor[1770] = 0.99935;
1063    TPCScaleFactor[1773] = 0.998531;
1064    TPCScaleFactor[1786] = 0.999125;
1065    TPCScaleFactor[1787] = 0.998479;
1066    TPCScaleFactor[1801] = 0.9775;
1067    TPCScaleFactor[1802] = 0.999095;
1068    TPCScaleFactor[1811] = 0.998197;
1069    TPCScaleFactor[1815] = 1.00413;
1070    TPCScaleFactor[1879] = 0.990916;
1071    TPCScaleFactor[1881] = 0.987241;
1072    TPCScaleFactor[2186] = 1.00048;
1073    TPCScaleFactor[2187] = 1.00057;
1074    TPCScaleFactor[2254] = 1.00588;
1075    TPCScaleFactor[2256] = 0.991503;
1076    TPCScaleFactor[2321] = 1.00057;
1077
1078
1079   // set all missing values to the value of the run before it ....
1080   for (int i=0; i<(fHighRunN-fLowRunN); i++) {    
1081     if (V0MScaleFactor[i] == 0.0) {     
1082       if (i==0) {
1083         V0MScaleFactor[i] = 1.00;
1084       } else {
1085         // search for last run number with non-zero value
1086         for (int j=i-1;j>=0;j--) {
1087           if (V0MScaleFactor[j] != 0.0) {
1088             V0MScaleFactor[i] = V0MScaleFactor[j];
1089             break;
1090           }
1091         }
1092       }
1093     }
1094   } // end loop over checking all run-numbers 
1095
1096   for (int i=0; i<(fHighRunN-fLowRunN); i++) {    
1097     if (SPDScaleFactor[i] == 0.0) {     
1098       if (i==0) {
1099         SPDScaleFactor[i] = 1.00;
1100       } else {
1101         for (int j=i-1;j>=0;j--) {
1102           if (SPDScaleFactor[j] != 0.0) {
1103             SPDScaleFactor[i] = SPDScaleFactor[j];
1104             break;
1105           }
1106         }
1107       }
1108     }
1109   } 
1110
1111   for (int i=0; i<(fHighRunN-fLowRunN); i++) {    
1112     if (TPCScaleFactor[i] == 0.0) {     
1113       if (i==0) {
1114         TPCScaleFactor[i] = 1.00;
1115       } else {
1116         for (int j=i-1;j>=0;j--) {
1117           if (TPCScaleFactor[j] != 0.0) {
1118             TPCScaleFactor[i] = TPCScaleFactor[j];
1119             break;
1120           }
1121         }
1122       }
1123     }
1124   } 
1125   
1126
1127   //    for (int i=0; i<(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << V0MScaleFactor[i] << " for Run " << i+fLowRunN << endl;
1128   
1129   return;
1130
1131 }
1132
1133
1134 //________________________________________________________________________  
1135 void AliCentralitySelectionTask::MyInitScaleFactorMC() 
1136 {
1137 // Initialize the MC scaling factors
1138   for (int i=0; i<(fHighRunN-fLowRunN); i++) V0MScaleFactorMC[i] = 0.0;
1139   // scale factors determined from <V0 charge> on a run-by-run basis
1140   V0MScaleFactorMC[0] = 0.75108;
1141   // set all missing values to the value of the run before it ....
1142   for (int i=0; i<(fHighRunN-fLowRunN); i++) {    
1143     if (V0MScaleFactorMC[i] == 0.0) {     
1144       if (i==0) {
1145         V0MScaleFactorMC[i] = 1.00;
1146       } else {
1147         // search for last run number with non-zero value
1148         for (int j=i-1;j>=0;j--) {
1149           if (V0MScaleFactorMC[j] != 0.0) {
1150             V0MScaleFactorMC[i] = V0MScaleFactorMC[j];
1151             break;
1152           }
1153         }
1154       }
1155     }
1156   } // end loop over checking all run-numbers 
1157
1158
1159   //    for (int i=0; i<(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << V0MScaleFactorMC[i] << " for Run " << i+fLowRunN << endl;
1160   
1161   return;
1162
1163 }
1164