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