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