Restoring previous fixes that were lost during one of the last commits in this class...
[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 <TTree.h>
23 #include <TList.h>
24 #include <TH1F.h>
25 #include <TH2F.h>
26 #include <TProfile.h>
27 #include <TFile.h>
28 #include <TObjString.h>
29 #include <TString.h>
30 #include <TCanvas.h>
31 #include <TROOT.h>
32 #include <TDirectory.h>
33 #include <TSystem.h>
34 #include <iostream>
35
36 #include "AliAnalysisManager.h"
37 #include "AliVEvent.h"
38 #include "AliESD.h"
39 #include "AliESDEvent.h"
40 #include "AliESDHeader.h"
41 #include "AliESDInputHandler.h"
42 #include "AliESDZDC.h"
43 #include "AliESDFMD.h"
44 #include "AliESDVZERO.h"
45 #include "AliESDCentrality.h"
46 #include "AliESDtrackCuts.h"
47 #include "AliMultiplicity.h"
48 #include "AliAODHandler.h"
49 #include "AliAODEvent.h"
50 #include "AliESDVertex.h"
51 #include "AliAODVertex.h"
52 #include "AliAODMCHeader.h"
53 #include "AliMCEvent.h"
54 #include "AliMCEventHandler.h"
55 #include "AliMCParticle.h"
56 #include "AliStack.h"
57 #include "AliHeader.h"
58 #include "AliAODMCParticle.h"
59 #include "AliAnalysisTaskSE.h"
60 #include "AliGenEventHeader.h"
61 #include "AliGenHijingEventHeader.h"
62 #include "AliPhysicsSelectionTask.h"
63 #include "AliPhysicsSelection.h"
64 #include "AliBackgroundSelection.h"
65 #include "AliCentralitySelectionTask.h"
66
67 ClassImp(AliCentralitySelectionTask)
68
69
70 //________________________________________________________________________
71 AliCentralitySelectionTask::AliCentralitySelectionTask():
72 AliAnalysisTaskSE(),
73   fDebug(0),
74   fAnalysisInput("ESD"),
75   fIsMCInput(kFALSE),
76   fFile(0),
77   fFile2(0),
78   fCurrentRun(-1),
79   fRunNo(-1),
80   fTrackCuts(0),
81   fCentV0M(0),
82   fCentFMD(0),
83   fCentTRK(0),
84   fCentTKL(0),
85   fCentCL0(0),
86   fCentCL1(0),
87   fCentV0MvsFMD(0),
88   fCentTKLvsV0M(0),
89   fCentZEMvsZDC(0),
90   fHtempV0M(0),
91   fHtempFMD(0),
92   fHtempTRK(0),
93   fHtempTKL(0),
94   fHtempCL0(0),
95   fHtempCL1(0),
96   fHtempV0MvsFMD(0),
97   fHtempTKLvsV0M(0),
98   fHtempZEMvsZDC(0),
99   fOutputList(0),
100   fHOutCentV0M     (0),
101   fHOutCentFMD     (0),
102   fHOutCentTRK     (0),
103   fHOutCentTKL     (0),
104   fHOutCentCL0     (0),
105   fHOutCentCL1     (0),
106   fHOutCentV0MvsFMD(0),
107   fHOutCentTKLvsV0M(0),
108   fHOutCentZEMvsZDC(0),
109   fHOutMultV0M(0),
110   fHOutMultFMD(0),
111   fHOutMultTRK(0),
112   fHOutMultTKL(0),
113   fHOutMultCL0(0),
114   fHOutMultCL1(0),
115   fHOutMultV0MvsZDC(0),
116   fHOutMultZEMvsZDC(0),
117   fHOutMultV0MvsCL1(0),
118   fHOutMultV0MvsTRK(0),
119   fHOutMultTRKvsCL1(0)
120 {   
121   // Default constructor
122   AliInfo("Centrality Selection enabled.");
123 }   
124
125 //________________________________________________________________________
126 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
127   AliAnalysisTaskSE(name),
128   fDebug(0),
129   fAnalysisInput("ESD"),
130   fIsMCInput(kFALSE),
131   fFile(0),
132   fFile2(0),
133   fCurrentRun(-1),
134   fRunNo(-1),
135   fTrackCuts(0),
136   fCentV0M(0),
137   fCentFMD(0),
138   fCentTRK(0),
139   fCentTKL(0),
140   fCentCL0(0),
141   fCentCL1(0),
142   fCentV0MvsFMD(0),
143   fCentTKLvsV0M(0),
144   fCentZEMvsZDC(0),
145   fHtempV0M(0),
146   fHtempFMD(0),
147   fHtempTRK(0),
148   fHtempTKL(0),
149   fHtempCL0(0),
150   fHtempCL1(0),
151   fHtempV0MvsFMD(0),
152   fHtempTKLvsV0M(0),
153   fHtempZEMvsZDC(0),
154   fOutputList(0),
155   fHOutCentV0M     (0),
156   fHOutCentFMD     (0),
157   fHOutCentTRK     (0),
158   fHOutCentTKL     (0),
159   fHOutCentCL0     (0),
160   fHOutCentCL1     (0),
161   fHOutCentV0MvsFMD(0),
162   fHOutCentTKLvsV0M(0),
163   fHOutCentZEMvsZDC(0),
164   fHOutMultV0M(0),
165   fHOutMultFMD(0),
166   fHOutMultTRK(0),
167   fHOutMultTKL(0),
168   fHOutMultCL0(0),
169   fHOutMultCL1(0),
170   fHOutMultV0MvsZDC(0),
171   fHOutMultZEMvsZDC(0),
172   fHOutMultV0MvsCL1(0),
173   fHOutMultV0MvsTRK(0),
174   fHOutMultTRKvsCL1(0)
175 {
176   // Default constructor
177   AliInfo("Centrality Selection enabled.");
178   DefineOutput(1, TList::Class());
179 }
180
181 //________________________________________________________________________
182 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
183 {
184   // Assignment operator
185   if (this!=&c) {
186     AliAnalysisTaskSE::operator=(c);
187   }
188   return *this;
189 }
190
191 //________________________________________________________________________
192 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
193   AliAnalysisTaskSE(ana),
194   fDebug(ana.fDebug),     
195   fAnalysisInput(ana.fDebug),
196   fIsMCInput(ana.fIsMCInput),
197   fFile(ana.fFile),
198   fFile2(ana.fFile2),
199   fCurrentRun(ana.fCurrentRun),
200   fRunNo(ana.fRunNo),
201   fTrackCuts(ana.fTrackCuts),
202   fCentV0M(ana.fCentV0M),
203   fCentFMD(ana.fCentFMD),
204   fCentTRK(ana.fCentTRK),
205   fCentTKL(ana.fCentTKL),
206   fCentCL0(ana.fCentCL0),
207   fCentCL1(ana.fCentCL1),
208   fCentV0MvsFMD(ana.fCentV0MvsFMD),
209   fCentTKLvsV0M(ana.fCentTKLvsV0M),
210   fCentZEMvsZDC(ana.fCentZEMvsZDC),
211   fHtempV0M(ana.fHtempV0M),
212   fHtempFMD(ana.fHtempFMD),
213   fHtempTRK(ana.fHtempTRK),
214   fHtempTKL(ana.fHtempTKL),
215   fHtempCL0(ana.fHtempCL0),
216   fHtempCL1(ana.fHtempCL1),
217   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
218   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
219   fHtempZEMvsZDC(ana.fHtempZEMvsZDC),
220   fOutputList(ana.fOutputList),
221   fHOutCentV0M     (ana.fHOutCentV0M     ),
222   fHOutCentFMD     (ana.fHOutCentFMD     ),
223   fHOutCentTRK     (ana.fHOutCentTRK     ),
224   fHOutCentTKL     (ana.fHOutCentTKL     ),
225   fHOutCentCL0     (ana.fHOutCentCL0     ),
226   fHOutCentCL1     (ana.fHOutCentCL1     ),
227   fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD),
228   fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M),
229   fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC),
230   fHOutMultV0M(ana.fHOutMultV0M),
231   fHOutMultFMD(ana.fHOutMultFMD),
232   fHOutMultTRK(ana.fHOutMultTRK),
233   fHOutMultTKL(ana.fHOutMultTKL),
234   fHOutMultCL0(ana.fHOutMultCL0),
235   fHOutMultCL1(ana.fHOutMultCL1),
236   fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC),
237   fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),
238   fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),
239   fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),
240   fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1)
241 {
242   // Copy Constructor   
243 }
244  
245 //________________________________________________________________________
246 AliCentralitySelectionTask::~AliCentralitySelectionTask()
247 {
248   // Destructor  
249   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
250   if (fTrackCuts) delete fTrackCuts;
251 }  
252
253 //________________________________________________________________________
254 void AliCentralitySelectionTask::UserCreateOutputObjects()
255 {  
256   // Create the output containers
257   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
258   AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
259
260   fOutputList = new TList();
261   fOutputList->SetOwner();
262   fHOutCentV0M     = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",101,-0.5,100.5);
263   fHOutCentFMD     = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",101,-0.5,100.5);
264   fHOutCentTRK     = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",101,-0.5,100.5);
265   fHOutCentTKL     = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",101,-0.5,100.5);
266   fHOutCentCL0     = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",101,-0.5,100.5);
267   fHOutCentCL1     = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",101,-0.5,100.5);
268   fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",101,-0.5,100.5);
269   fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",101,-0.5,100.5);
270   fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",101,-0.5,100.5);
271
272   fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,25000);
273   fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000);
274   fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000);
275   fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000);
276   fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000);
277   fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000);
278   fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,25000,500,0,6000);
279   fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,6000);
280   fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
281   fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
282   fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
283
284   fOutputList->Add(  fHOutCentV0M     );
285   fOutputList->Add(  fHOutCentFMD     );
286   fOutputList->Add(  fHOutCentTRK     );
287   fOutputList->Add(  fHOutCentTKL     );
288   fOutputList->Add(  fHOutCentCL0     );
289   fOutputList->Add(  fHOutCentCL1     );
290   fOutputList->Add(  fHOutCentV0MvsFMD);
291   fOutputList->Add(  fHOutCentTKLvsV0M);
292   fOutputList->Add(  fHOutCentZEMvsZDC);
293   fOutputList->Add(  fHOutMultV0M); 
294   fOutputList->Add(  fHOutMultFMD); 
295   fOutputList->Add(  fHOutMultTRK); 
296   fOutputList->Add(  fHOutMultTKL); 
297   fOutputList->Add(  fHOutMultCL0); 
298   fOutputList->Add(  fHOutMultCL1); 
299   fOutputList->Add(  fHOutMultV0MvsZDC);
300   fOutputList->Add(  fHOutMultZEMvsZDC);
301   fOutputList->Add(  fHOutMultV0MvsCL1);
302   fOutputList->Add(  fHOutMultV0MvsTRK);
303   fOutputList->Add(  fHOutMultTRKvsCL1);
304
305
306   fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
307
308   PostData(1, fOutputList); 
309 }
310
311 //________________________________________________________________________
312 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
313
314   // Execute analysis for current event:
315   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
316   
317   Float_t  zncEnergy = 0.;          //  ZNC Energy
318   Float_t  zpcEnergy = 0.;          //  ZPC Energy
319   Float_t  znaEnergy = 0.;          //  ZNA Energy
320   Float_t  zpaEnergy = 0.;          //  ZPA Energy
321   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
322   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
323   
324   Int_t    nTracks = 0;             //  no. tracks
325   Int_t    nTracklets = 0;          //  no. tracklets
326   Int_t    nClusters[6];            //  no. clusters on 6 ITS layers
327   Int_t    nChips[2];               //  no. chips on 2 SPD layers
328   Float_t  spdCorr =0;              //  corrected spd2 multiplicity
329
330   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
331   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
332   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
333   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
334
335   Short_t v0Corr = 0;               // corrected V0 multiplicity
336   Short_t v0CorrResc = 0;           // corrected and rescaled V0 multiplicity
337
338   Float_t zvtx =0;                  // z-vertex SPD
339  
340   AliESDCentrality *esdCent = 0;
341
342   if(fAnalysisInput.CompareTo("ESD")==0){
343
344     AliVEvent* event = InputEvent();
345     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
346   
347     if (fRunNo<=0) {
348       if (SetupRun(esd)<0)
349          AliFatal("Centrality File not available for this run");
350     }
351
352     esdCent = esd->GetCentrality();
353
354     // ***** V0 info    
355     AliESDVZERO* esdV0 = esd->GetVZEROData();
356     multV0A=esdV0->GetMTotV0A();
357     multV0C=esdV0->GetMTotV0C();
358
359     float v0CorrR;
360     v0Corr = (Short_t)GetCorrV0(esd,v0CorrR);
361     v0CorrResc = (Short_t)v0CorrR;
362
363     // ***** Vertex Info
364     const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
365     zvtx        = vtxESD->GetZ(); 
366
367     // ***** CB info (tracklets, clusters, chips)
368     //nTracks    = event->GetNumberOfTracks();     
369     nTracks    = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
370
371     const AliMultiplicity *mult = esd->GetMultiplicity();
372
373     nTracklets = mult->GetNumberOfTracklets();
374
375     for(Int_t ilay=0; ilay<6; ilay++){
376       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
377     }
378
379     for(Int_t ilay=0; ilay<2; ilay++){
380       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
381     }
382
383     spdCorr = GetCorrSPD2(nClusters[1],zvtx);    
384
385     // ***** FMD info
386     AliESDFMD *fmd = esd->GetFMDData();
387     Float_t totalMultA = 0;
388     Float_t totalMultC = 0;
389     const Float_t fFMDLowCut = 0.4;
390     
391     for(UShort_t det=1;det<=3;det++) {
392       Int_t nRings = (det==1 ? 1 : 2);
393       for (UShort_t ir = 0; ir < nRings; ir++) {          
394         Char_t   ring = (ir == 0 ? 'I' : 'O');
395         UShort_t nsec = (ir == 0 ? 20  : 40);
396         UShort_t nstr = (ir == 0 ? 512 : 256);
397         for(UShort_t sec =0; sec < nsec;  sec++)  {
398           for(UShort_t strip = 0; strip < nstr; strip++) {
399             
400             Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
401             if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
402             
403             Float_t nParticles=0;
404             
405             if(FMDmult > fFMDLowCut) {
406               nParticles = 1.;
407             }
408             
409             if (det<3) totalMultA = totalMultA + nParticles;
410             else totalMultC = totalMultC + nParticles;
411             
412           }
413         }
414       }
415     }
416     multFMDA = totalMultA;
417     multFMDC = totalMultC;
418     
419     // ***** ZDC info
420     AliESDZDC *esdZDC = esd->GetESDZDC();
421     zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
422     zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
423     znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
424     zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
425     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
426     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
427     
428   }   
429   else if(fAnalysisInput.CompareTo("AOD")==0){
430     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
431     // to be implemented
432     printf("  AOD analysis not yet implemented!!!\n\n");
433     return;
434   }
435
436   // ***** Centrality Selection
437   if(fHtempV0M)  fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
438   ///  else     printf("  Centrality by V0 not available!!!\n\n");
439   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
440   //  else     printf("  Centrality by FMD not available!!!\n\n");
441   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
442   //  else     printf("  Centrality by TRK not available!!!\n\n");
443   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
444   //  else     printf("  Centrality by TKL not available!!!\n\n");
445   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
446   //  else     printf("  Centrality by CL0 not available!!!\n\n");
447   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
448   ///  else     printf("  Centrality by CL1 not available!!!\n\n");
449   
450   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
451   //  else     printf("  Centrality by V0 vs FMD not available!!!\n\n");
452   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
453   //  else     printf("  Centrality by V0 vs TKL not available!!!\n\n");
454   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
455   //  else     printf("  Centrality by ZEM vs ZDC not available!!!\n\n");
456
457   esdCent->SetCentralityV0M(fCentV0M);
458   esdCent->SetCentralityFMD(fCentFMD);
459   esdCent->SetCentralityTRK(fCentTRK);
460   esdCent->SetCentralityTKL(fCentTKL);
461   esdCent->SetCentralityCL0(fCentCL0);
462   esdCent->SetCentralityCL1(fCentCL1);
463   esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
464   esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
465   esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
466
467   fHOutCentV0M->Fill(fCentV0M);
468   fHOutCentFMD->Fill(fCentFMD);
469   fHOutCentTRK->Fill(fCentTRK);
470   fHOutCentTKL->Fill(fCentTKL);
471   fHOutCentCL0->Fill(fCentCL0);
472   fHOutCentCL1->Fill(fCentCL1);
473   fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
474   fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
475   fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
476   fHOutMultV0M->Fill(v0Corr);
477   fHOutMultFMD->Fill((multFMDA+multFMDC));
478   fHOutMultTRK->Fill(nTracks);
479   fHOutMultTKL->Fill(nTracklets);
480   fHOutMultCL0->Fill(nClusters[0]);
481   fHOutMultCL1->Fill(spdCorr);
482   fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
483   fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
484   fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
485   fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
486   fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
487
488   PostData(1, fOutputList); 
489 }
490
491 //________________________________________________________________________
492 void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename) 
493 {
494   //  Read centrality histograms
495   TDirectory *owd = gDirectory;
496   // Check if the file is present
497   TString path = gSystem->ExpandPathName(fCentfilename.Data());
498   if (gSystem->AccessPathName(path)) {
499      AliError(Form("File %s does not exist", path.Data()));
500      return;
501   }
502   fFile  = TFile::Open(fCentfilename);
503   owd->cd();
504   fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));
505   fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));
506   fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));
507   fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));
508   fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));
509   fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));
510   owd->cd();
511 }  
512
513 //________________________________________________________________________
514 void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2) 
515 {
516   //  Read centrality histograms
517   TDirectory *owd = gDirectory;
518   TString path = gSystem->ExpandPathName(fCentfilename2.Data());
519   if (gSystem->AccessPathName(path)) {
520      AliError(Form("File %s does not exist", path.Data()));
521      return;
522   }   
523   fFile2  = TFile::Open(fCentfilename2);
524   owd->cd();
525   fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
526   fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
527   fHtempZEMvsZDC  = (TH1F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
528   owd->cd();
529 }
530
531 //________________________________________________________________________
532 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
533 {
534   // Terminate analysis
535   if (fFile && fFile->IsOpen())
536     fFile->Close();  
537   if (fFile2 && fFile2->IsOpen())
538     fFile2->Close();  
539 }
540 //________________________________________________________________________
541 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* esd)
542 {
543   // Setup files for run
544
545   if (!esd)
546     return -1;
547
548   // check if something to be done
549   if (fCurrentRun == esd->GetRunNumber())
550     return 0;
551   else
552     fCurrentRun = esd->GetRunNumber();
553   
554   AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
555
556   fRunNo = fCurrentRun;
557
558   // CHANGE HERE FOR RUN RANGES
559   if ( fRunNo <= 137162 ) fRunNo = 137161;
560   else if ( fRunNo == 137365 ) fRunNo = 137366;
561   else if ( fRunNo >= 137366 ) fRunNo = 137366;
562   // CHANGE HERE FOR RUN RANGES
563   
564   TString fileName(Form("$ALICE_ROOT/ANALYSIS/macros/AliCentralityBy1D_%d.root", fRunNo));
565   TString fileName2(Form("$ALICE_ROOT/ANALYSIS/macros/AliCentralityByFunction_%d.root", fRunNo));
566   
567   AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));
568   ReadCentralityHistos(fileName.Data());
569   ReadCentralityHistos2(fileName2.Data());
570   if (!fFile && !fFile2) {
571      AliFatal(Form("Run %d not known to centrality selection!", fCurrentRun));       
572      return -1;
573   }   
574   return 0;
575 }
576 //________________________________________________________________________
577 Float_t AliCentralitySelectionTask::GetCorrV0(const AliESDEvent* esd, float &v0CorrResc) 
578 {
579   // correct V0 non-linearity, prepare a version rescaled to SPD2 corr
580   Double_t *par0;
581   Double_t *par1;
582   Double_t *par2;
583   
584   Double_t par0_137161[64] = { 6.71e-02 , 6.86e-02 , 7.06e-02 , 6.32e-02 , 
585                         5.91e-02 , 6.07e-02 , 5.78e-02 , 5.73e-02 , 5.91e-02 , 6.22e-02 , 
586                         5.90e-02 , 6.11e-02 , 5.55e-02 , 5.29e-02 , 5.19e-02 , 5.56e-02 , 
587                         6.25e-02 , 7.03e-02 , 5.64e-02 , 5.81e-02 , 4.57e-02 , 5.30e-02 , 
588                         5.13e-02 , 6.43e-02 , 6.27e-02 , 6.48e-02 , 6.07e-02 , 1.01e-01 , 
589                         6.68e-02 , 7.16e-02 , 6.36e-02 , 5.95e-02 , 2.52e-02 , 2.82e-02 , 
590                         2.56e-02 , 2.86e-02 , 2.82e-02 , 2.10e-02 , 2.13e-02 , 2.32e-02 , 
591                         2.75e-02 , 4.34e-02 , 3.78e-02 , 4.52e-02 , 4.11e-02 , 3.89e-02 , 
592                         4.10e-02 , 3.73e-02 , 4.51e-02 , 5.07e-02 , 5.42e-02 , 4.74e-02 , 
593                         4.33e-02 , 4.44e-02 , 4.64e-02 , 3.01e-02 , 6.38e-02 , 5.26e-02 , 
594                         4.99e-02 , 5.26e-02 , 5.47e-02 , 3.84e-02 , 5.00e-02 , 5.20e-02 };
595   Double_t par1_137161[64] = { -6.68e-05 , -7.78e-05 , -6.88e-05 , -5.92e-05 , 
596                         -2.43e-05 , -3.54e-05 , -2.91e-05 , -1.99e-05 , -1.40e-05 , -4.01e-05 , 
597                         -2.29e-05 , -3.68e-05 , -2.53e-05 , -2.44e-06 , -9.22e-06 , -1.51e-05 , 
598                         -2.80e-05 , -2.34e-05 , -1.72e-05 , -1.81e-05 , -1.29e-05 , -2.65e-05 , 
599                         -1.61e-05 , -2.86e-05 , -1.74e-05 , -4.23e-05 , -3.41e-05 , -1.05e-04 , 
600                         -2.76e-05 , -4.71e-05 , -3.06e-05 , -2.32e-05 , -1.55e-06 , 2.15e-05 , 
601                         1.40e-05 , 2.16e-05 , 1.21e-05 , 3.05e-06 , 1.67e-05 , -3.84e-06 , 
602                         3.09e-06 , 1.50e-05 , 3.47e-06 , 4.87e-06 , -3.71e-07 , -1.75e-06 , 
603                         -1.80e-06 , 9.99e-06 , -6.46e-06 , -4.91e-06 , 1.33e-05 , -2.52e-07 , 
604                         -3.85e-06 , 4.94e-06 , -2.48e-07 , -1.20e-05 , 2.07e-06 , 6.12e-06 , 
605                         -1.18e-06 , 4.54e-06 , -1.54e-05 , -1.25e-05 , 1.46e-06 , -6.67e-06 };
606   Double_t par2_137161[64] = { 1.29e-08 , 1.51e-08 , 1.43e-08 , 1.11e-08 , 
607                         5.04e-09 , 6.99e-09 , 5.58e-09 , 4.15e-09 , 4.00e-09 , 8.22e-09 , 
608                         4.97e-09 , 7.66e-09 , 4.91e-09 , 1.10e-09 , 2.64e-09 , 3.64e-09 , 
609                         5.76e-09 , 5.46e-09 , 3.38e-09 , 3.47e-09 , 2.43e-09 , 4.13e-09 , 
610                         2.80e-09 , 5.80e-09 , 3.86e-09 , 7.46e-09 , 5.98e-09 , 2.58e-08 , 
611                         5.50e-09 , 8.72e-09 , 5.23e-09 , 4.37e-09 , 2.33e-09 , -6.01e-10 , 
612                         3.99e-11 , -2.02e-10 , 7.67e-10 , 2.03e-09 , 1.17e-10 , 2.56e-09 , 
613                         1.16e-09 , -4.75e-10 , 1.28e-09 , 1.23e-09 , 1.62e-09 , 1.61e-09 , 
614                         1.93e-09 , 2.97e-10 , 2.21e-09 , 2.16e-09 , 5.22e-10 , 1.03e-09 , 
615                         1.56e-09 , 5.00e-10 , 1.01e-09 , 2.93e-09 , 1.05e-09 , 9.96e-11 , 
616                         1.21e-09 , 7.45e-10 , 3.07e-09 , 2.31e-09 , 6.70e-10 , 1.89e-09 };
617
618   Double_t par0_137366[64] = { 7.12e-02 , 7.34e-02 , 7.39e-02 , 6.54e-02 , 6.11e-02 , 6.31e-02 , 6.15e-02 , 
619                                6.00e-02 , 6.10e-02 , 6.49e-02 , 6.17e-02 , 6.33e-02 , 6.00e-02 , 5.48e-02 , 
620                                5.44e-02 , 5.81e-02 , 6.49e-02 , 7.07e-02 , 5.91e-02 , 6.18e-02 , 4.82e-02 , 
621                                5.67e-02 , 5.36e-02 , 6.60e-02 , 6.37e-02 , 6.78e-02 , 6.31e-02 , 1.04e-01 , 
622                                6.91e-02 , 7.32e-02 , 6.61e-02 , 6.16e-02 , 2.64e-02 , 2.81e-02 , 2.64e-02 , 
623                                2.85e-02 , 2.87e-02 , 2.18e-02 , 2.19e-02 , 2.43e-02 , 2.81e-02 , 4.37e-02 , 
624                                3.90e-02 , 4.66e-02 , 4.24e-02 , 4.09e-02 , 4.21e-02 , 3.88e-02 , 4.83e-02 , 
625                                5.23e-02 , 5.44e-02 , 4.85e-02 , 4.42e-02 , 4.58e-02 , 4.74e-02 , 3.14e-02 , 
626                                6.31e-02 , 5.30e-02 , 5.01e-02 , 5.33e-02 , 5.70e-02 , 3.95e-02 , 4.98e-02 , 5.31e-02 };
627   Double_t par1_137366[64] = { -6.99e-05 , -6.99e-05 , -6.94e-05 , -6.55e-05 , -3.55e-05 , -4.50e-05 , 
628                                -3.10e-05 , -2.81e-05 , -2.29e-05 , -3.89e-05 , -2.53e-05 , -4.25e-05 ,
629                                -1.87e-05 , -2.01e-05 , -1.53e-05 , -2.14e-05 , -2.86e-05 , -4.70e-05 ,
630                                -2.23e-05 , -3.30e-05 ,-9.74e-06 , -2.62e-05 , -1.76e-05 , -2.38e-05 , 
631                                -2.40e-05 , -3.43e-05 , -2.75e-05 , -6.86e-05 ,-2.35e-05 , -4.45e-05 , 
632                                -2.51e-05 , -2.20e-05 , -1.25e-16 , -2.04e-17 , -2.06e-17 , -3.74e-19 ,
633                                -1.18e-18 , -2.02e-15 , -3.78e-06 , -1.26e-06 , -2.71e-06 , -6.23e-17 , 
634                                -7.39e-08 , -1.76e-16 , -8.98e-06 , -4.10e-18 , -1.34e-05 , -1.06e-16 , 
635                                -3.34e-06 , -1.04e-05 , -5.28e-06 , -7.34e-06 , -1.05e-05 , -7.68e-06 ,
636                                -1.78e-05 , -1.19e-05 , -1.78e-05 , -1.34e-06 , -9.23e-06 , -3.34e-06 ,
637                                -8.02e-06 , -1.39e-05 , -1.38e-05 , -1.40e-05 };
638   Double_t par2_137366[64] = { 1.41e-08 , 1.47e-08 , 1.48e-08 , 1.24e-08 , 6.82e-09 , 8.73e-09 , 6.26e-09 , 
639                                5.53e-09 , 5.40e-09 , 7.93e-09 , 5.49e-09 , 8.77e-09 , 4.21e-09 , 3.93e-09 , 
640                                3.60e-09 , 4.67e-09 , 5.59e-09 , 8.81e-09 , 3.89e-09 , 6.19e-09 , 1.97e-09 , 
641                                4.38e-09 , 3.26e-09 , 5.00e-09 , 4.58e-09 , 6.39e-09 , 5.03e-09 , 1.30e-08 , 
642                                4.95e-09 , 8.26e-09 , 4.57e-09 , 4.10e-09 , 2.35e-09 , 2.30e-09 , 2.15e-09 , 
643                                2.27e-09 , 2.17e-09 , 2.27e-09 , 2.97e-09 , 2.25e-09 , 1.69e-09 , 1.44e-09 , 
644                                1.66e-09 , 1.75e-09 , 2.88e-09 , 1.82e-09 , 3.64e-09 , 1.80e-09 , 1.71e-09 , 
645                                2.66e-09 , 3.01e-09 , 1.95e-09 , 2.64e-09 , 2.42e-09 , 3.68e-09 , 2.66e-09 , 
646                                3.92e-09 , 1.18e-09 , 2.26e-09 , 1.57e-09 , 2.02e-09 , 2.71e-09 , 2.99e-09 , 3.04e-09 }; 
647   
648   
649   if (esd->GetRunNumber() <= 137165) {
650     par0=par0_137161;
651     par1=par1_137161;
652     par2=par2_137161;
653   }  else  {
654     par0=par0_137366;
655     par1=par1_137366;
656     par2=par2_137366;
657  }
658   //
659   Float_t multCorr = 0;
660   Float_t multCorr2 = 0;
661   Float_t multChCorr[64];
662   AliESDVZERO* esdV0 = esd->GetVZEROData();
663   for(Int_t i = 0; i < 64; ++i) {
664     Double_t b = (esdV0->GetMultiplicity(i)*par1[i]-par0[i]);
665     Double_t s = (b*b-4.*par2[i]*esdV0->GetMultiplicity(i)*esdV0->GetMultiplicity(i));
666     Double_t n = (s<0) ? -b : (-b + TMath::Sqrt(s));
667     multChCorr[i] = 2.*esdV0->GetMultiplicity(i)/n*par0[i];
668     multCorr += multChCorr[i];
669     multCorr2 += (multChCorr[i]/par0[i]/64.);
670   }
671   v0CorrResc =  multCorr2;
672   return multCorr;
673 }
674
675 //____________________________________________________________________
676 Float_t AliCentralitySelectionTask::GetCorrSPD2(Float_t spd2raw,Float_t zv)
677 {
678   // renormalize N spd2 clusters at given Zv to acceptance at Zv=0
679   const double pars[] = {8.10030e-01,-2.80364e-03,-7.19504e-04};
680   zv -= pars[0];
681   float corr = 1 + zv*(pars[1] + zv*pars[2]);
682   return corr>0 ? spd2raw/corr : -1;
683 }