5b06b7e9cc3e20d9c3560e7307d301be3f60fe3c
[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 = GetCorrV0(esd,v0CorrR,fRunNo);
361     v0Corr = (Short_t)v0Corr;
362     v0CorrResc = (Short_t)v0CorrR;
363
364     // ***** Vertex Info
365     const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
366     zvtx        = vtxESD->GetZ(); 
367
368     // ***** CB info (tracklets, clusters, chips)
369     //nTracks    = event->GetNumberOfTracks();     
370     nTracks    = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
371
372     const AliMultiplicity *mult = esd->GetMultiplicity();
373
374     nTracklets = mult->GetNumberOfTracklets();
375
376     for(Int_t ilay=0; ilay<6; ilay++){
377       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
378     }
379
380     for(Int_t ilay=0; ilay<2; ilay++){
381       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
382     }
383
384     spdCorr = GetCorrSPD2(nClusters[1],zvtx);    
385
386     // ***** FMD info
387     AliESDFMD *fmd = esd->GetFMDData();
388     Float_t totalMultA = 0;
389     Float_t totalMultC = 0;
390     const Float_t fFMDLowCut = 0.4;
391     
392     for(UShort_t det=1;det<=3;det++) {
393       Int_t nRings = (det==1 ? 1 : 2);
394       for (UShort_t ir = 0; ir < nRings; ir++) {          
395         Char_t   ring = (ir == 0 ? 'I' : 'O');
396         UShort_t nsec = (ir == 0 ? 20  : 40);
397         UShort_t nstr = (ir == 0 ? 512 : 256);
398         for(UShort_t sec =0; sec < nsec;  sec++)  {
399           for(UShort_t strip = 0; strip < nstr; strip++) {
400             
401             Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
402             if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
403             
404             Float_t nParticles=0;
405             
406             if(FMDmult > fFMDLowCut) {
407               nParticles = 1.;
408             }
409             
410             if (det<3) totalMultA = totalMultA + nParticles;
411             else totalMultC = totalMultC + nParticles;
412             
413           }
414         }
415       }
416     }
417     multFMDA = totalMultA;
418     multFMDC = totalMultC;
419     
420     // ***** ZDC info
421     AliESDZDC *esdZDC = esd->GetESDZDC();
422     zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
423     zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
424     znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
425     zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
426     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
427     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
428     
429   }   
430   else if(fAnalysisInput.CompareTo("AOD")==0){
431     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
432     // to be implemented
433     printf("  AOD analysis not yet implemented!!!\n\n");
434     return;
435   }
436
437   // ***** Centrality Selection
438   if(fHtempV0M)  fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
439   ///  else     printf("  Centrality by V0 not available!!!\n\n");
440   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
441   //  else     printf("  Centrality by FMD not available!!!\n\n");
442   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
443   //  else     printf("  Centrality by TRK not available!!!\n\n");
444   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
445   //  else     printf("  Centrality by TKL not available!!!\n\n");
446   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
447   //  else     printf("  Centrality by CL0 not available!!!\n\n");
448   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
449   ///  else     printf("  Centrality by CL1 not available!!!\n\n");
450   
451   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
452   //  else     printf("  Centrality by V0 vs FMD not available!!!\n\n");
453   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
454   //  else     printf("  Centrality by V0 vs TKL not available!!!\n\n");
455   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
456   //  else     printf("  Centrality by ZEM vs ZDC not available!!!\n\n");
457
458   esdCent->SetCentralityV0M(fCentV0M);
459   esdCent->SetCentralityFMD(fCentFMD);
460   esdCent->SetCentralityTRK(fCentTRK);
461   esdCent->SetCentralityTKL(fCentTKL);
462   esdCent->SetCentralityCL0(fCentCL0);
463   esdCent->SetCentralityCL1(fCentCL1);
464   esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
465   esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
466   esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
467
468   fHOutCentV0M->Fill(fCentV0M);
469   fHOutCentFMD->Fill(fCentFMD);
470   fHOutCentTRK->Fill(fCentTRK);
471   fHOutCentTKL->Fill(fCentTKL);
472   fHOutCentCL0->Fill(fCentCL0);
473   fHOutCentCL1->Fill(fCentCL1);
474   fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
475   fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
476   fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
477   fHOutMultV0M->Fill(v0Corr);
478   fHOutMultFMD->Fill((multFMDA+multFMDC));
479   fHOutMultTRK->Fill(nTracks);
480   fHOutMultTKL->Fill(nTracklets);
481   fHOutMultCL0->Fill(nClusters[0]);
482   fHOutMultCL1->Fill(spdCorr);
483   fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
484   fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
485   fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
486   fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
487   fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
488
489   PostData(1, fOutputList); 
490 }
491
492 //________________________________________________________________________
493 void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename) 
494 {
495   //  Read centrality histograms
496   TDirectory *owd = gDirectory;
497   // Check if the file is present
498   TString path = gSystem->ExpandPathName(fCentfilename.Data());
499   if (gSystem->AccessPathName(path)) {
500      AliError(Form("File %s does not exist", path.Data()));
501      return;
502   }
503   fFile  = TFile::Open(fCentfilename);
504   owd->cd();
505   fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));
506   fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));
507   fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));
508   fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));
509   fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));
510   fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));
511   owd->cd();
512 }  
513
514 //________________________________________________________________________
515 void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2) 
516 {
517   //  Read centrality histograms
518   TDirectory *owd = gDirectory;
519   TString path = gSystem->ExpandPathName(fCentfilename2.Data());
520   if (gSystem->AccessPathName(path)) {
521      AliError(Form("File %s does not exist", path.Data()));
522      return;
523   }   
524   fFile2  = TFile::Open(fCentfilename2);
525   owd->cd();
526   fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
527   fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
528   fHtempZEMvsZDC  = (TH1F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
529   owd->cd();
530 }
531
532 //________________________________________________________________________
533 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
534 {
535   // Terminate analysis
536   if (fFile && fFile->IsOpen())
537     fFile->Close();  
538   if (fFile2 && fFile2->IsOpen())
539     fFile2->Close();  
540 }
541 //________________________________________________________________________
542 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* esd)
543 {
544   // Setup files for run
545
546   if (!esd)
547     return -1;
548
549   // check if something to be done
550   if (fCurrentRun == esd->GetRunNumber())
551     return 0;
552   else
553     fCurrentRun = esd->GetRunNumber();
554   
555   AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
556
557   fRunNo = fCurrentRun;
558
559   // CHANGE HERE FOR RUN RANGES
560   if ( fRunNo <= 137162 ) fRunNo = 137161;
561   else if ( fRunNo == 137365 ) fRunNo = 137366;
562   else if ( fRunNo >= 137366 ) fRunNo = 137366;
563   // CHANGE HERE FOR RUN RANGES
564   
565   TString fileName(Form("$ALICE_ROOT/ANALYSIS/macros/AliCentralityBy1D_%d.root", fRunNo));
566   TString fileName2(Form("$ALICE_ROOT/ANALYSIS/macros/AliCentralityByFunction_%d.root", fRunNo));
567   
568   AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));
569   ReadCentralityHistos(fileName.Data());
570   ReadCentralityHistos2(fileName2.Data());
571   if (!fFile && !fFile2) {
572      AliFatal(Form("Run %d not known to centrality selection!", fCurrentRun));       
573      return -1;
574   }   
575   return 0;
576 }
577 //________________________________________________________________________
578 Float_t AliCentralitySelectionTask::GetCorrV0(const AliESDEvent* esd, float &v0CorrResc, int run) 
579 {
580   // correct V0 non-linearity, prepare a version rescaled to SPD2 corr
581   Double_t *par0;
582   Double_t *par1;
583   Double_t *par2;
584   
585   Double_t par0_137161[64] = { 6.71e-02 , 6.86e-02 , 7.06e-02 , 6.32e-02 , 
586                         5.91e-02 , 6.07e-02 , 5.78e-02 , 5.73e-02 , 5.91e-02 , 6.22e-02 , 
587                         5.90e-02 , 6.11e-02 , 5.55e-02 , 5.29e-02 , 5.19e-02 , 5.56e-02 , 
588                         6.25e-02 , 7.03e-02 , 5.64e-02 , 5.81e-02 , 4.57e-02 , 5.30e-02 , 
589                         5.13e-02 , 6.43e-02 , 6.27e-02 , 6.48e-02 , 6.07e-02 , 1.01e-01 , 
590                         6.68e-02 , 7.16e-02 , 6.36e-02 , 5.95e-02 , 2.52e-02 , 2.82e-02 , 
591                         2.56e-02 , 2.86e-02 , 2.82e-02 , 2.10e-02 , 2.13e-02 , 2.32e-02 , 
592                         2.75e-02 , 4.34e-02 , 3.78e-02 , 4.52e-02 , 4.11e-02 , 3.89e-02 , 
593                         4.10e-02 , 3.73e-02 , 4.51e-02 , 5.07e-02 , 5.42e-02 , 4.74e-02 , 
594                         4.33e-02 , 4.44e-02 , 4.64e-02 , 3.01e-02 , 6.38e-02 , 5.26e-02 , 
595                         4.99e-02 , 5.26e-02 , 5.47e-02 , 3.84e-02 , 5.00e-02 , 5.20e-02 };
596   Double_t par1_137161[64] = { -6.68e-05 , -7.78e-05 , -6.88e-05 , -5.92e-05 , 
597                         -2.43e-05 , -3.54e-05 , -2.91e-05 , -1.99e-05 , -1.40e-05 , -4.01e-05 , 
598                         -2.29e-05 , -3.68e-05 , -2.53e-05 , -2.44e-06 , -9.22e-06 , -1.51e-05 , 
599                         -2.80e-05 , -2.34e-05 , -1.72e-05 , -1.81e-05 , -1.29e-05 , -2.65e-05 , 
600                         -1.61e-05 , -2.86e-05 , -1.74e-05 , -4.23e-05 , -3.41e-05 , -1.05e-04 , 
601                         -2.76e-05 , -4.71e-05 , -3.06e-05 , -2.32e-05 , -1.55e-06 , 2.15e-05 , 
602                         1.40e-05 , 2.16e-05 , 1.21e-05 , 3.05e-06 , 1.67e-05 , -3.84e-06 , 
603                         3.09e-06 , 1.50e-05 , 3.47e-06 , 4.87e-06 , -3.71e-07 , -1.75e-06 , 
604                         -1.80e-06 , 9.99e-06 , -6.46e-06 , -4.91e-06 , 1.33e-05 , -2.52e-07 , 
605                         -3.85e-06 , 4.94e-06 , -2.48e-07 , -1.20e-05 , 2.07e-06 , 6.12e-06 , 
606                         -1.18e-06 , 4.54e-06 , -1.54e-05 , -1.25e-05 , 1.46e-06 , -6.67e-06 };
607   Double_t par2_137161[64] = { 1.29e-08 , 1.51e-08 , 1.43e-08 , 1.11e-08 , 
608                         5.04e-09 , 6.99e-09 , 5.58e-09 , 4.15e-09 , 4.00e-09 , 8.22e-09 , 
609                         4.97e-09 , 7.66e-09 , 4.91e-09 , 1.10e-09 , 2.64e-09 , 3.64e-09 , 
610                         5.76e-09 , 5.46e-09 , 3.38e-09 , 3.47e-09 , 2.43e-09 , 4.13e-09 , 
611                         2.80e-09 , 5.80e-09 , 3.86e-09 , 7.46e-09 , 5.98e-09 , 2.58e-08 , 
612                         5.50e-09 , 8.72e-09 , 5.23e-09 , 4.37e-09 , 2.33e-09 , -6.01e-10 , 
613                         3.99e-11 , -2.02e-10 , 7.67e-10 , 2.03e-09 , 1.17e-10 , 2.56e-09 , 
614                         1.16e-09 , -4.75e-10 , 1.28e-09 , 1.23e-09 , 1.62e-09 , 1.61e-09 , 
615                         1.93e-09 , 2.97e-10 , 2.21e-09 , 2.16e-09 , 5.22e-10 , 1.03e-09 , 
616                         1.56e-09 , 5.00e-10 , 1.01e-09 , 2.93e-09 , 1.05e-09 , 9.96e-11 , 
617                         1.21e-09 , 7.45e-10 , 3.07e-09 , 2.31e-09 , 6.70e-10 , 1.89e-09 };
618
619   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 , 
620                                6.00e-02 , 6.10e-02 , 6.49e-02 , 6.17e-02 , 6.33e-02 , 6.00e-02 , 5.48e-02 , 
621                                5.44e-02 , 5.81e-02 , 6.49e-02 , 7.07e-02 , 5.91e-02 , 6.18e-02 , 4.82e-02 , 
622                                5.67e-02 , 5.36e-02 , 6.60e-02 , 6.37e-02 , 6.78e-02 , 6.31e-02 , 1.04e-01 , 
623                                6.91e-02 , 7.32e-02 , 6.61e-02 , 6.16e-02 , 2.64e-02 , 2.81e-02 , 2.64e-02 , 
624                                2.85e-02 , 2.87e-02 , 2.18e-02 , 2.19e-02 , 2.43e-02 , 2.81e-02 , 4.37e-02 , 
625                                3.90e-02 , 4.66e-02 , 4.24e-02 , 4.09e-02 , 4.21e-02 , 3.88e-02 , 4.83e-02 , 
626                                5.23e-02 , 5.44e-02 , 4.85e-02 , 4.42e-02 , 4.58e-02 , 4.74e-02 , 3.14e-02 , 
627                                6.31e-02 , 5.30e-02 , 5.01e-02 , 5.33e-02 , 5.70e-02 , 3.95e-02 , 4.98e-02 , 5.31e-02 };
628   Double_t par1_137366[64] = { -6.99e-05 , -6.99e-05 , -6.94e-05 , -6.55e-05 , -3.55e-05 , -4.50e-05 , 
629                                -3.10e-05 , -2.81e-05 , -2.29e-05 , -3.89e-05 , -2.53e-05 , -4.25e-05 ,
630                                -1.87e-05 , -2.01e-05 , -1.53e-05 , -2.14e-05 , -2.86e-05 , -4.70e-05 ,
631                                -2.23e-05 , -3.30e-05 ,-9.74e-06 , -2.62e-05 , -1.76e-05 , -2.38e-05 , 
632                                -2.40e-05 , -3.43e-05 , -2.75e-05 , -6.86e-05 ,-2.35e-05 , -4.45e-05 , 
633                                -2.51e-05 , -2.20e-05 , -1.25e-16 , -2.04e-17 , -2.06e-17 , -3.74e-19 ,
634                                -1.18e-18 , -2.02e-15 , -3.78e-06 , -1.26e-06 , -2.71e-06 , -6.23e-17 , 
635                                -7.39e-08 , -1.76e-16 , -8.98e-06 , -4.10e-18 , -1.34e-05 , -1.06e-16 , 
636                                -3.34e-06 , -1.04e-05 , -5.28e-06 , -7.34e-06 , -1.05e-05 , -7.68e-06 ,
637                                -1.78e-05 , -1.19e-05 , -1.78e-05 , -1.34e-06 , -9.23e-06 , -3.34e-06 ,
638                                -8.02e-06 , -1.39e-05 , -1.38e-05 , -1.40e-05 };
639   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 , 
640                                5.53e-09 , 5.40e-09 , 7.93e-09 , 5.49e-09 , 8.77e-09 , 4.21e-09 , 3.93e-09 , 
641                                3.60e-09 , 4.67e-09 , 5.59e-09 , 8.81e-09 , 3.89e-09 , 6.19e-09 , 1.97e-09 , 
642                                4.38e-09 , 3.26e-09 , 5.00e-09 , 4.58e-09 , 6.39e-09 , 5.03e-09 , 1.30e-08 , 
643                                4.95e-09 , 8.26e-09 , 4.57e-09 , 4.10e-09 , 2.35e-09 , 2.30e-09 , 2.15e-09 , 
644                                2.27e-09 , 2.17e-09 , 2.27e-09 , 2.97e-09 , 2.25e-09 , 1.69e-09 , 1.44e-09 , 
645                                1.66e-09 , 1.75e-09 , 2.88e-09 , 1.82e-09 , 3.64e-09 , 1.80e-09 , 1.71e-09 , 
646                                2.66e-09 , 3.01e-09 , 1.95e-09 , 2.64e-09 , 2.42e-09 , 3.68e-09 , 2.66e-09 , 
647                                3.92e-09 , 1.18e-09 , 2.26e-09 , 1.57e-09 , 2.02e-09 , 2.71e-09 , 2.99e-09 , 3.04e-09 }; 
648   
649   
650   if (run==137161) {
651     par0=par0_137161;
652     par1=par1_137161;
653     par2=par2_137161;
654   }  else  {
655     par0=par0_137366;
656     par1=par1_137366;
657     par2=par2_137366;
658  }
659   //
660   Float_t multCorr = 0;
661   Float_t multCorr2 = 0;
662   Float_t multChCorr[64];
663   AliESDVZERO* esdV0 = esd->GetVZEROData();
664   for(Int_t i = 0; i < 64; ++i) {
665     Double_t b = (esdV0->GetMultiplicity(i)*par1[i]-par0[i]);
666     Double_t s = (b*b-4.*par2[i]*esdV0->GetMultiplicity(i)*esdV0->GetMultiplicity(i));
667     Double_t n;
668     if (s<0) {
669       printf("FPE %d %.2f %.2f %.2e\n",i,esdV0->GetMultiplicity(i),b,(b*b-4.*par2[i]*esdV0->GetMultiplicity(i)*esdV0->GetMultiplicity(i)));
670       n = -b;
671     }
672     else {
673       n = (-b + TMath::Sqrt(s));
674     }
675     multChCorr[i] = 2.*esdV0->GetMultiplicity(i)/n*par0[i];
676     multCorr += multChCorr[i];
677     multCorr2 += (multChCorr[i]/par0[i]/64.);
678   }
679   v0CorrResc =  multCorr2;
680   return multCorr;
681 }
682
683 //____________________________________________________________________
684 Float_t AliCentralitySelectionTask::GetCorrSPD2(Float_t spd2raw,Float_t zv)
685 {
686   // renormalize N spd2 clusters at given Zv to acceptance at Zv=0
687   const double pars[] = {8.10030e-01,-2.80364e-03,-7.19504e-04};
688   zv -= pars[0];
689   float corr = 1 + zv*(pars[1] + zv*pars[2]);
690   return corr>0 ? spd2raw/corr : -1;
691 }