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