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