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