git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Restoring previous fixes that were lost during one of the last commits in this class...
[u/mrichter/AliRoot.git] / ANALYSIS / AliCentralitySelectionTask.cxx
CommitLineData
e8472cd0 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>
68e344b7 28#include <TObjString.h>
e8472cd0 29#include <TString.h>
30#include <TCanvas.h>
9f14d90a 31#include <TROOT.h>
32#include <TDirectory.h>
4fd15e3c 33#include <TSystem.h>
e8472cd0 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"
05915818 46#include "AliESDtrackCuts.h"
e8472cd0 47#include "AliMultiplicity.h"
48#include "AliAODHandler.h"
49#include "AliAODEvent.h"
62164ba5 50#include "AliESDVertex.h"
e8472cd0 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
67ClassImp(AliCentralitySelectionTask)
68
69
70//________________________________________________________________________
71AliCentralitySelectionTask::AliCentralitySelectionTask():
72AliAnalysisTaskSE(),
73 fDebug(0),
74 fAnalysisInput("ESD"),
75 fIsMCInput(kFALSE),
76 fFile(0),
d15bf53f 77 fFile2(0),
68e344b7 78 fCurrentRun(-1),
62164ba5 79 fRunNo(-1),
05915818 80 fTrackCuts(0),
d15bf53f 81 fCentV0M(0),
82 fCentFMD(0),
83 fCentTRK(0),
84 fCentTKL(0),
85 fCentCL0(0),
be0d4e9b 86 fCentCL1(0),
d15bf53f 87 fCentV0MvsFMD(0),
88 fCentTKLvsV0M(0),
89 fCentZEMvsZDC(0),
90 fHtempV0M(0),
91 fHtempFMD(0),
92 fHtempTRK(0),
93 fHtempTKL(0),
94 fHtempCL0(0),
be0d4e9b 95 fHtempCL1(0),
d15bf53f 96 fHtempV0MvsFMD(0),
97 fHtempTKLvsV0M(0),
05915818 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),
36ee40df 108 fHOutCentZEMvsZDC(0),
109 fHOutMultV0M(0),
110 fHOutMultFMD(0),
111 fHOutMultTRK(0),
112 fHOutMultTKL(0),
113 fHOutMultCL0(0),
114 fHOutMultCL1(0),
115 fHOutMultV0MvsZDC(0),
ab57f513 116 fHOutMultZEMvsZDC(0),
117 fHOutMultV0MvsCL1(0),
118 fHOutMultV0MvsTRK(0),
119 fHOutMultTRKvsCL1(0)
e8472cd0 120{
121 // Default constructor
68e344b7 122 AliInfo("Centrality Selection enabled.");
e8472cd0 123}
124
125//________________________________________________________________________
126AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
127 AliAnalysisTaskSE(name),
128 fDebug(0),
129 fAnalysisInput("ESD"),
130 fIsMCInput(kFALSE),
131 fFile(0),
d15bf53f 132 fFile2(0),
68e344b7 133 fCurrentRun(-1),
62164ba5 134 fRunNo(-1),
05915818 135 fTrackCuts(0),
d15bf53f 136 fCentV0M(0),
137 fCentFMD(0),
138 fCentTRK(0),
139 fCentTKL(0),
140 fCentCL0(0),
be0d4e9b 141 fCentCL1(0),
d15bf53f 142 fCentV0MvsFMD(0),
143 fCentTKLvsV0M(0),
144 fCentZEMvsZDC(0),
145 fHtempV0M(0),
146 fHtempFMD(0),
147 fHtempTRK(0),
148 fHtempTKL(0),
149 fHtempCL0(0),
be0d4e9b 150 fHtempCL1(0),
d15bf53f 151 fHtempV0MvsFMD(0),
152 fHtempTKLvsV0M(0),
05915818 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),
36ee40df 163 fHOutCentZEMvsZDC(0),
164 fHOutMultV0M(0),
165 fHOutMultFMD(0),
166 fHOutMultTRK(0),
167 fHOutMultTKL(0),
168 fHOutMultCL0(0),
169 fHOutMultCL1(0),
170 fHOutMultV0MvsZDC(0),
ab57f513 171 fHOutMultZEMvsZDC(0),
172 fHOutMultV0MvsCL1(0),
173 fHOutMultV0MvsTRK(0),
174 fHOutMultTRKvsCL1(0)
e8472cd0 175{
176 // Default constructor
68e344b7 177 AliInfo("Centrality Selection enabled.");
05915818 178 DefineOutput(1, TList::Class());
e8472cd0 179}
180
181//________________________________________________________________________
182AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
183{
184 // Assignment operator
185 if (this!=&c) {
186 AliAnalysisTaskSE::operator=(c);
187 }
188 return *this;
189}
190
191//________________________________________________________________________
192AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
193 AliAnalysisTaskSE(ana),
194 fDebug(ana.fDebug),
195 fAnalysisInput(ana.fDebug),
196 fIsMCInput(ana.fIsMCInput),
197 fFile(ana.fFile),
d15bf53f 198 fFile2(ana.fFile2),
68e344b7 199 fCurrentRun(ana.fCurrentRun),
62164ba5 200 fRunNo(ana.fRunNo),
05915818 201 fTrackCuts(ana.fTrackCuts),
d15bf53f 202 fCentV0M(ana.fCentV0M),
203 fCentFMD(ana.fCentFMD),
204 fCentTRK(ana.fCentTRK),
205 fCentTKL(ana.fCentTKL),
206 fCentCL0(ana.fCentCL0),
be0d4e9b 207 fCentCL1(ana.fCentCL1),
d15bf53f 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),
be0d4e9b 216 fHtempCL1(ana.fHtempCL1),
d15bf53f 217 fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
218 fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
05915818 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),
36ee40df 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),
ab57f513 237 fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),
238 fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),
239 fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),
240 fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1)
e8472cd0 241{
242 // Copy Constructor
243}
244
245//________________________________________________________________________
68e344b7 246AliCentralitySelectionTask::~AliCentralitySelectionTask()
247{
05915818 248 // Destructor
249 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
4fd15e3c 250 if (fTrackCuts) delete fTrackCuts;
68e344b7 251}
e8472cd0 252
253//________________________________________________________________________
254void AliCentralitySelectionTask::UserCreateOutputObjects()
255{
256 // Create the output containers
257 if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
68e344b7 258 AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
05915818 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
36ee40df 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);
ab57f513 278 fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,25000,500,0,6000);
a7d618ed 279 fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,6000);
ab57f513 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);
36ee40df 283
05915818 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);
36ee40df 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);
ab57f513 301 fOutputList->Add( fHOutMultV0MvsCL1);
302 fOutputList->Add( fHOutMultV0MvsTRK);
303 fOutputList->Add( fHOutMultTRKvsCL1);
36ee40df 304
305
4fd15e3c 306 fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
05915818 307
308 PostData(1, fOutputList);
e8472cd0 309}
310
311//________________________________________________________________________
312void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
313{
314 // Execute analysis for current event:
315 if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
316
36ee40df 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
d15bf53f 321 Float_t zem1Energy = 0.; // ZEM1 Energy
322 Float_t zem2Energy = 0.; // ZEM2 Energy
e8472cd0 323
be0d4e9b 324 Int_t nTracks = 0; // no. tracks
d15bf53f 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
62164ba5 328 Float_t spdCorr =0; // corrected spd2 multiplicity
e8472cd0 329
d15bf53f 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
e8472cd0 334
62164ba5 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
e8472cd0 340 AliESDCentrality *esdCent = 0;
341
342 if(fAnalysisInput.CompareTo("ESD")==0){
343
344 AliVEvent* event = InputEvent();
345 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
05915818 346
347 if (fRunNo<=0) {
348 if (SetupRun(esd)<0)
349 AliFatal("Centrality File not available for this run");
350 }
68e344b7 351
e8472cd0 352 esdCent = esd->GetCentrality();
353
354 // ***** V0 info
355 AliESDVZERO* esdV0 = esd->GetVZEROData();
356 multV0A=esdV0->GetMTotV0A();
357 multV0C=esdV0->GetMTotV0C();
62164ba5 358
359 float v0CorrR;
0caa916b 360 v0Corr = (Short_t)GetCorrV0(esd,v0CorrR);
62164ba5 361 v0CorrResc = (Short_t)v0CorrR;
362
363 // ***** Vertex Info
364 const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
365 zvtx = vtxESD->GetZ();
366
e8472cd0 367 // ***** CB info (tracklets, clusters, chips)
05915818 368 //nTracks = event->GetNumberOfTracks();
05915818 369 nTracks = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
e8472cd0 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 }
62164ba5 382
383 spdCorr = GetCorrSPD2(nClusters[1],zvtx);
e8472cd0 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();
36ee40df 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.;
e8472cd0 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
d15bf53f 436 // ***** Centrality Selection
62164ba5 437 if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
68e344b7 438 /// else printf(" Centrality by V0 not available!!!\n\n");
99df2fab 439 if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
68e344b7 440 // else printf(" Centrality by FMD not available!!!\n\n");
99df2fab 441 if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
68e344b7 442 // else printf(" Centrality by TRK not available!!!\n\n");
99df2fab 443 if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
68e344b7 444 // else printf(" Centrality by TKL not available!!!\n\n");
99df2fab 445 if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
68e344b7 446 // else printf(" Centrality by CL0 not available!!!\n\n");
62164ba5 447 if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
68e344b7 448 /// else printf(" Centrality by CL1 not available!!!\n\n");
d15bf53f 449
99df2fab 450 if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
68e344b7 451 // else printf(" Centrality by V0 vs FMD not available!!!\n\n");
99df2fab 452 if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
68e344b7 453 // else printf(" Centrality by V0 vs TKL not available!!!\n\n");
99df2fab 454 if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
68e344b7 455 // else printf(" Centrality by ZEM vs ZDC not available!!!\n\n");
456
d15bf53f 457 esdCent->SetCentralityV0M(fCentV0M);
458 esdCent->SetCentralityFMD(fCentFMD);
459 esdCent->SetCentralityTRK(fCentTRK);
460 esdCent->SetCentralityTKL(fCentTKL);
461 esdCent->SetCentralityCL0(fCentCL0);
be0d4e9b 462 esdCent->SetCentralityCL1(fCentCL1);
d15bf53f 463 esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
464 esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
465 esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
68e344b7 466
05915818 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);
36ee40df 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));
ab57f513 484 fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
485 fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
486 fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
05915818 487
488 PostData(1, fOutputList);
d15bf53f 489}
490
68e344b7 491//________________________________________________________________________
05915818 492void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename)
d15bf53f 493{
05915818 494 // Read centrality histograms
495 TDirectory *owd = gDirectory;
4fd15e3c 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 }
05915818 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"));
4fd15e3c 510 owd->cd();
05915818 511}
d15bf53f 512
68e344b7 513//________________________________________________________________________
05915818 514void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2)
d15bf53f 515{
05915818 516 // Read centrality histograms
517 TDirectory *owd = gDirectory;
4fd15e3c 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 }
05915818 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"));
4fd15e3c 528 owd->cd();
e8472cd0 529}
530
531//________________________________________________________________________
532void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
533{
534 // Terminate analysis
68e344b7 535 if (fFile && fFile->IsOpen())
536 fFile->Close();
537 if (fFile2 && fFile2->IsOpen())
538 fFile2->Close();
e8472cd0 539}
540//________________________________________________________________________
68e344b7 541Int_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
62164ba5 556 fRunNo = fCurrentRun;
68e344b7 557
558 // CHANGE HERE FOR RUN RANGES
05915818 559 if ( fRunNo <= 137162 ) fRunNo = 137161;
62164ba5 560 else if ( fRunNo == 137365 ) fRunNo = 137366;
05915818 561 else if ( fRunNo >= 137366 ) fRunNo = 137366;
68e344b7 562 // CHANGE HERE FOR RUN RANGES
05915818 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));
05915818 566
4fd15e3c 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;
68e344b7 575}
62164ba5 576//________________________________________________________________________
0caa916b 577Float_t AliCentralitySelectionTask::GetCorrV0(const AliESDEvent* esd, float &v0CorrResc)
62164ba5 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
0caa916b 649 if (esd->GetRunNumber() <= 137165) {
62164ba5 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));
0caa916b 666 Double_t n = (s<0) ? -b : (-b + TMath::Sqrt(s));
62164ba5 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//____________________________________________________________________
676Float_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}