1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //*****************************************************
17 // Class AliCentralitySelectionTask
18 // Class to analyze determine centrality
19 // author: Alberica Toia
20 //*****************************************************
22 #include "AliCentralitySelectionTask.h"
31 #include <TObjString.h>
35 #include <TDirectory.h>
39 #include "AliAnalysisManager.h"
40 #include "AliHeader.h"
41 #include "AliVEvent.h"
43 #include "AliESDEvent.h"
44 #include "AliESDHeader.h"
45 #include "AliESDInputHandler.h"
46 #include "AliESDZDC.h"
47 #include "AliESDFMD.h"
48 #include "AliESDVZERO.h"
49 #include "AliESDtrackCuts.h"
50 #include "AliESDVertex.h"
51 #include "AliCentrality.h"
52 #include "AliOADBCentrality.h"
53 #include "AliOADBContainer.h"
54 #include "AliMultiplicity.h"
55 #include "AliAODHandler.h"
56 #include "AliAODHeader.h"
57 #include "AliAODEvent.h"
58 #include "AliAODVertex.h"
59 #include "AliAODVZERO.h"
60 #include "AliAODTracklets.h"
61 #include "AliAODMCHeader.h"
62 #include "AliMCEventHandler.h"
63 #include "AliMCEvent.h"
64 #include "AliAODMCParticle.h"
65 #include "AliMCParticle.h"
67 #include "AliAnalysisTaskSE.h"
68 #include "AliGenEventHeader.h"
69 #include "AliGenHijingEventHeader.h"
70 #include "AliPhysicsSelectionTask.h"
71 #include "AliPhysicsSelection.h"
72 #include "AliBackgroundSelection.h"
73 #include "AliESDUtils.h"
75 ClassImp(AliCentralitySelectionTask)
78 //________________________________________________________________________
79 AliCentralitySelectionTask::AliCentralitySelectionTask():
81 fAnalysisInput("ESD"),
91 fV0MSPDOutlierPar0(0),
92 fV0MSPDOutlierPar1(0),
93 fV0MTPCOutlierPar0(0),
94 fV0MTPCOutlierPar1(0),
95 fV0MSPDSigmaOutlierPar0(0),
96 fV0MSPDSigmaOutlierPar1(0),
97 fV0MSPDSigmaOutlierPar2(0),
98 fV0MTPCSigmaOutlierPar0(0),
99 fV0MTPCSigmaOutlierPar1(0),
100 fV0MTPCSigmaOutlierPar2(0),
101 fV0MZDCOutlierPar0(0),
102 fV0MZDCOutlierPar1(0),
103 fV0MZDCEcalOutlierPar0(0),
104 fV0MZDCEcalOutlierPar1(0),
134 fHOutCentV0MvsFMD(0),
135 fHOutCentTKLvsV0M(0),
136 fHOutCentZEMvsZDC(0),
137 fHOutCentV0MvsCentCL1(0),
138 fHOutCentV0MvsCentTRK(0),
139 fHOutCentTRKvsCentCL1(0),
140 fHOutCentV0MvsCentZDC(0),
148 fHOutMultV0MvsZDN(0),
149 fHOutMultZEMvsZDN(0),
150 fHOutMultV0MvsZDC(0),
151 fHOutMultZEMvsZDC(0),
152 fHOutMultZEMvsZDCw(0),
153 fHOutMultV0MvsCL1(0),
154 fHOutMultV0MvsTRK(0),
155 fHOutMultTRKvsCL1(0),
156 fHOutMultV0MvsV0O(0),
157 fHOutMultV0OvsCL1(0),
158 fHOutMultV0OvsTRK(0),
159 fHOutCentV0Mqual1(0),
160 fHOutCentTRKqual1(0),
161 fHOutCentCL1qual1(0),
162 fHOutMultV0MvsCL1qual1(0),
163 fHOutMultV0MvsTRKqual1(0),
164 fHOutMultTRKvsCL1qual1(0),
165 fHOutCentV0Mqual2(0),
166 fHOutCentTRKqual2(0),
167 fHOutCentCL1qual2(0),
168 fHOutMultV0MvsCL1qual2(0),
169 fHOutMultV0MvsTRKqual2(0),
170 fHOutMultTRKvsCL1qual2(0),
174 // Default constructor
175 AliInfo("Centrality Selection enabled.");
179 fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO."
180 ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
183 //________________________________________________________________________
184 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
185 AliAnalysisTaskSE(name),
186 fAnalysisInput("ESD"),
195 fV0MScaleFactorMC(0),
196 fV0MSPDOutlierPar0(0),
197 fV0MSPDOutlierPar1(0),
198 fV0MTPCOutlierPar0(0),
199 fV0MTPCOutlierPar1(0),
200 fV0MSPDSigmaOutlierPar0(0),
201 fV0MSPDSigmaOutlierPar1(0),
202 fV0MSPDSigmaOutlierPar2(0),
203 fV0MTPCSigmaOutlierPar0(0),
204 fV0MTPCSigmaOutlierPar1(0),
205 fV0MTPCSigmaOutlierPar2(0),
206 fV0MZDCOutlierPar0(0),
207 fV0MZDCOutlierPar1(0),
208 fV0MZDCEcalOutlierPar0(0),
209 fV0MZDCEcalOutlierPar1(0),
239 fHOutCentV0MvsFMD(0),
240 fHOutCentTKLvsV0M(0),
241 fHOutCentZEMvsZDC(0),
242 fHOutCentV0MvsCentCL1(0),
243 fHOutCentV0MvsCentTRK(0),
244 fHOutCentTRKvsCentCL1(0),
245 fHOutCentV0MvsCentZDC(0),
253 fHOutMultV0MvsZDN(0),
254 fHOutMultZEMvsZDN(0),
255 fHOutMultV0MvsZDC(0),
256 fHOutMultZEMvsZDC(0),
257 fHOutMultZEMvsZDCw(0),
258 fHOutMultV0MvsCL1(0),
259 fHOutMultV0MvsTRK(0),
260 fHOutMultTRKvsCL1(0),
261 fHOutMultV0MvsV0O(0),
262 fHOutMultV0OvsCL1(0),
263 fHOutMultV0OvsTRK(0),
264 fHOutCentV0Mqual1(0),
265 fHOutCentTRKqual1(0),
266 fHOutCentCL1qual1(0),
267 fHOutMultV0MvsCL1qual1(0),
268 fHOutMultV0MvsTRKqual1(0),
269 fHOutMultTRKvsCL1qual1(0),
270 fHOutCentV0Mqual2(0),
271 fHOutCentTRKqual2(0),
272 fHOutCentCL1qual2(0),
273 fHOutMultV0MvsCL1qual2(0),
274 fHOutMultV0MvsTRKqual2(0),
275 fHOutMultTRKvsCL1qual2(0),
279 // Default constructor
280 AliInfo("Centrality Selection enabled.");
281 DefineOutput(1, TList::Class());
284 fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO."
285 ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
288 //________________________________________________________________________
289 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
291 // Assignment operator
293 AliAnalysisTaskSE::operator=(c);
298 //________________________________________________________________________
299 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
300 AliAnalysisTaskSE(ana),
301 fAnalysisInput(ana.fDebug),
302 fIsMCInput(ana.fIsMCInput),
304 fCurrentRun(ana.fCurrentRun),
305 fUseScaling(ana.fUseScaling),
306 fUseCleaning(ana.fUseCleaning),
307 fV0MScaleFactor(ana.fV0MScaleFactor),
308 fSPDScaleFactor(ana.fSPDScaleFactor),
309 fTPCScaleFactor(ana.fTPCScaleFactor),
310 fV0MScaleFactorMC(ana.fV0MScaleFactorMC),
311 fV0MSPDOutlierPar0(ana.fV0MSPDOutlierPar0),
312 fV0MSPDOutlierPar1(ana.fV0MSPDOutlierPar1),
313 fV0MTPCOutlierPar0(ana.fV0MTPCOutlierPar0),
314 fV0MTPCOutlierPar1(ana.fV0MTPCOutlierPar1),
315 fV0MSPDSigmaOutlierPar0(ana.fV0MSPDSigmaOutlierPar0),
316 fV0MSPDSigmaOutlierPar1(ana.fV0MSPDSigmaOutlierPar1),
317 fV0MSPDSigmaOutlierPar2(ana.fV0MSPDSigmaOutlierPar2),
318 fV0MTPCSigmaOutlierPar0(ana.fV0MTPCSigmaOutlierPar0),
319 fV0MTPCSigmaOutlierPar1(ana.fV0MTPCSigmaOutlierPar1),
320 fV0MTPCSigmaOutlierPar2(ana.fV0MTPCSigmaOutlierPar2),
321 fV0MZDCOutlierPar0(ana.fV0MZDCOutlierPar0),
322 fV0MZDCOutlierPar1(ana.fV0MZDCOutlierPar1),
323 fV0MZDCEcalOutlierPar0(ana.fV0MZDCEcalOutlierPar0),
324 fV0MZDCEcalOutlierPar1(ana.fV0MZDCEcalOutlierPar1),
325 fTrackCuts(ana.fTrackCuts),
327 fOutliersCut(ana.fOutliersCut),
328 fQuality(ana.fQuality),
329 fCentV0M(ana.fCentV0M),
330 fCentFMD(ana.fCentFMD),
331 fCentTRK(ana.fCentTRK),
332 fCentTKL(ana.fCentTKL),
333 fCentCL0(ana.fCentCL0),
334 fCentCL1(ana.fCentCL1),
335 fCentV0MvsFMD(ana.fCentV0MvsFMD),
336 fCentTKLvsV0M(ana.fCentTKLvsV0M),
337 fCentZEMvsZDC(ana.fCentZEMvsZDC),
338 fHtempV0M(ana.fHtempV0M),
339 fHtempFMD(ana.fHtempFMD),
340 fHtempTRK(ana.fHtempTRK),
341 fHtempTKL(ana.fHtempTKL),
342 fHtempCL0(ana.fHtempCL0),
343 fHtempCL1(ana.fHtempCL1),
344 fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
345 fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
346 fHtempZEMvsZDC(ana.fHtempZEMvsZDC),
347 fOutputList(ana.fOutputList),
348 fHOutCentV0M (ana.fHOutCentV0M ),
349 fHOutCentFMD (ana.fHOutCentFMD ),
350 fHOutCentTRK (ana.fHOutCentTRK ),
351 fHOutCentTKL (ana.fHOutCentTKL ),
352 fHOutCentCL0 (ana.fHOutCentCL0 ),
353 fHOutCentCL1 (ana.fHOutCentCL1 ),
354 fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD),
355 fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M),
356 fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC),
357 fHOutCentV0MvsCentCL1(ana.fHOutCentV0MvsCentCL1),
358 fHOutCentV0MvsCentTRK(ana.fHOutCentV0MvsCentTRK),
359 fHOutCentTRKvsCentCL1(ana.fHOutCentTRKvsCentCL1),
360 fHOutCentV0MvsCentZDC(ana.fHOutCentV0MvsCentZDC),
361 fHOutMultV0M(ana.fHOutMultV0M),
362 fHOutMultV0O(ana.fHOutMultV0O),
363 fHOutMultFMD(ana.fHOutMultFMD),
364 fHOutMultTRK(ana.fHOutMultTRK),
365 fHOutMultTKL(ana.fHOutMultTKL),
366 fHOutMultCL0(ana.fHOutMultCL0),
367 fHOutMultCL1(ana.fHOutMultCL1),
368 fHOutMultV0MvsZDN(ana.fHOutMultV0MvsZDN),
369 fHOutMultZEMvsZDN(ana.fHOutMultZEMvsZDN),
370 fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC),
371 fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),
372 fHOutMultZEMvsZDCw(ana.fHOutMultZEMvsZDCw),
373 fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),
374 fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),
375 fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1),
376 fHOutMultV0MvsV0O(ana.fHOutMultV0MvsV0O),
377 fHOutMultV0OvsCL1(ana.fHOutMultV0OvsCL1),
378 fHOutMultV0OvsTRK(ana.fHOutMultV0OvsTRK),
379 fHOutCentV0Mqual1(ana.fHOutCentV0Mqual1),
380 fHOutCentTRKqual1(ana.fHOutCentTRKqual1),
381 fHOutCentCL1qual1(ana.fHOutCentCL1qual1),
382 fHOutMultV0MvsCL1qual1(ana.fHOutMultV0MvsCL1qual1),
383 fHOutMultV0MvsTRKqual1(ana.fHOutMultV0MvsTRKqual1),
384 fHOutMultTRKvsCL1qual1(ana.fHOutMultTRKvsCL1qual1),
385 fHOutCentV0Mqual2(ana.fHOutCentV0Mqual2),
386 fHOutCentTRKqual2(ana.fHOutCentTRKqual2),
387 fHOutCentCL1qual2(ana.fHOutCentCL1qual2),
388 fHOutMultV0MvsCL1qual2(ana.fHOutMultV0MvsCL1qual2),
389 fHOutMultV0MvsTRKqual2(ana.fHOutMultV0MvsTRKqual2),
390 fHOutMultTRKvsCL1qual2(ana.fHOutMultTRKvsCL1qual2),
391 fHOutQuality(ana.fHOutQuality),
392 fHOutVertex(ana.fHOutVertex)
398 //________________________________________________________________________
399 AliCentralitySelectionTask::~AliCentralitySelectionTask()
402 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
403 if (fTrackCuts) delete fTrackCuts;
406 //________________________________________________________________________
407 void AliCentralitySelectionTask::UserCreateOutputObjects()
409 // Create the output containers
410 if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
411 AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
413 fOutputList = new TList();
414 fOutputList->SetOwner();
415 fHOutCentV0M = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",505,0,101);
416 fHOutCentFMD = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",505,0,101);
417 fHOutCentTRK = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",505,0,101);
418 fHOutCentTKL = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",505,0,101);
419 fHOutCentCL0 = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",505,0,101);
420 fHOutCentCL1 = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",505,0,101);
421 fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",505,0,101);
422 fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",505,0,101);
423 fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",505,0,101);
424 fHOutCentV0MvsCentCL1= new TH2F("fHOutCentV0MvsCentCL1","fHOutCentV0MvsCentCL1; Cent V0; Cent SPD",505,0,101,505,0,101);
425 fHOutCentV0MvsCentTRK= new TH2F("fHOutCentV0MvsCentTRK","fHOutCentV0MvsCentTRK; Cent V0; Cent TPC",505,0,101,505,0,101);
426 fHOutCentTRKvsCentCL1= new TH2F("fHOutCentTRKvsCentCL1","fHOutCentTRKvsCentCL1; Cent TPC; Cent SPD",505,0,101,505,0,101);
427 fHOutCentV0MvsCentZDC= new TH2F("fHOutCentV0MvsCentZDC","fHOutCentV0MvsCentZDC; Cent V0; Cent ZDC",505,0,101,505,0,101);
429 fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,30000);
430 fHOutMultV0O = new TH1F("fHOutMultV0O","fHOutMultV0O; Multiplicity V0",30000,0,30000);
431 fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000);
432 fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000);
433 fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000);
434 fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000);
435 fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000);
436 fHOutMultV0MvsZDN = new TH2F("fHOutMultV0MvsZDN","fHOutMultV0MvsZDN; Multiplicity V0; Energy ZDC-N",500,0,30000,500,0,180000);
437 fHOutMultZEMvsZDN = new TH2F("fHOutMultZEMvsZDN","fHOutMultZEMvsZDN; Energy ZEM; Energy ZDC-N",500,0,2500,500,0,180000);
438 fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,30000,500,0,200000);
439 fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,200000);
440 fHOutMultZEMvsZDCw = new TH2F("fHOutMultZEMvsZDCw","fHOutMultZEMvsZDCw; Energy ZEM; Energy ZDC (weigthed with V0 percentile)",500,0,2500,500,0,200000);
441 fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,30000,700,0,7000);
442 fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,30000,400,0,4000);
443 fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
444 fHOutMultV0MvsV0O = new TH2F("fHOutMultV0MvsV0O","fHOutMultV0MvsV0O; Multiplicity V0; Multiplicity V0 Online",500,0,30000,500,0,30000);
445 fHOutMultV0OvsCL1 = new TH2F("fHOutMultV0OvsCL1","fHOutMultV0OvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,30000,700,0,7000);
446 fHOutMultV0OvsTRK = new TH2F("fHOutMultV0OvsTRK","fHOutMultV0OvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,30000,400,0,4000);
448 fHOutCentV0Mqual1 = new TH1F("fHOutCentV0M_qual1","fHOutCentV0M_qual1; Centrality V0",505,0,101);
449 fHOutCentTRKqual1 = new TH1F("fHOutCentTRK_qual1","fHOutCentTRK_qual1; Centrality TPC",505,0,101);
450 fHOutCentCL1qual1 = new TH1F("fHOutCentCL1_qual1","fHOutCentCL1_qual1; Centrality SPD outer",505,0,101);
451 fHOutMultV0MvsCL1qual1 = new TH2F("fHOutMultV0MvsCL1_qual1","fHOutMultV0MvsCL1_qual1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
452 fHOutMultV0MvsTRKqual1 = new TH2F("fHOutMultV0MvsTRK_qual1","fHOutMultV0MvsTRK_qual1; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
453 fHOutMultTRKvsCL1qual1 = new TH2F("fHOutMultTRKvsCL1_qual1","fHOutMultTRKvsCL1_qual1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
455 fHOutCentV0Mqual2 = new TH1F("fHOutCentV0M_qual2","fHOutCentV0M_qual2; Centrality V0",505,0,101);
456 fHOutCentTRKqual2 = new TH1F("fHOutCentTRK_qual2","fHOutCentTRK_qual2; Centrality TPC",505,0,101);
457 fHOutCentCL1qual2 = new TH1F("fHOutCentCL1_qual2","fHOutCentCL1_qual2; Centrality SPD outer",505,0,101);
458 fHOutMultV0MvsCL1qual2 = new TH2F("fHOutMultV0MvsCL1_qual2","fHOutMultV0MvsCL1_qual2; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
459 fHOutMultV0MvsTRKqual2 = new TH2F("fHOutMultV0MvsTRK_qual2","fHOutMultV0MvsTRK_qual2; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
460 fHOutMultTRKvsCL1qual2 = new TH2F("fHOutMultTRKvsCL1_qual2","fHOutMultTRKvsCL1_qual2; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
462 fHOutQuality = new TH1F("fHOutQuality", "fHOutQuality", 100,-0.5,99.5);
463 fHOutVertex = new TH1F("fHOutVertex", "fHOutVertex", 100,-20,20);
465 fOutputList->Add( fHOutCentV0M );
466 fOutputList->Add( fHOutCentFMD );
467 fOutputList->Add( fHOutCentTRK );
468 fOutputList->Add( fHOutCentTKL );
469 fOutputList->Add( fHOutCentCL0 );
470 fOutputList->Add( fHOutCentCL1 );
471 fOutputList->Add( fHOutCentV0MvsFMD);
472 fOutputList->Add( fHOutCentTKLvsV0M);
473 fOutputList->Add( fHOutCentZEMvsZDC);
474 fOutputList->Add( fHOutCentV0MvsCentCL1);
475 fOutputList->Add( fHOutCentV0MvsCentTRK);
476 fOutputList->Add( fHOutCentTRKvsCentCL1);
477 fOutputList->Add( fHOutCentV0MvsCentZDC);
478 fOutputList->Add( fHOutMultV0M);
479 fOutputList->Add( fHOutMultV0O);
480 fOutputList->Add( fHOutMultFMD);
481 fOutputList->Add( fHOutMultTRK);
482 fOutputList->Add( fHOutMultTKL);
483 fOutputList->Add( fHOutMultCL0);
484 fOutputList->Add( fHOutMultCL1);
485 fOutputList->Add( fHOutMultV0MvsZDN);
486 fOutputList->Add( fHOutMultZEMvsZDN);
487 fOutputList->Add( fHOutMultV0MvsZDC);
488 fOutputList->Add( fHOutMultZEMvsZDC);
489 fOutputList->Add( fHOutMultZEMvsZDCw);
490 fOutputList->Add( fHOutMultV0MvsCL1);
491 fOutputList->Add( fHOutMultV0MvsTRK);
492 fOutputList->Add( fHOutMultTRKvsCL1);
493 fOutputList->Add( fHOutMultV0MvsV0O);
494 fOutputList->Add( fHOutMultV0OvsCL1);
495 fOutputList->Add( fHOutMultV0OvsTRK);
496 fOutputList->Add( fHOutCentV0Mqual1 );
497 fOutputList->Add( fHOutCentTRKqual1 );
498 fOutputList->Add( fHOutCentCL1qual1 );
499 fOutputList->Add( fHOutMultV0MvsCL1qual1);
500 fOutputList->Add( fHOutMultV0MvsTRKqual1);
501 fOutputList->Add( fHOutMultTRKvsCL1qual1);
502 fOutputList->Add( fHOutCentV0Mqual2 );
503 fOutputList->Add( fHOutCentTRKqual2 );
504 fOutputList->Add( fHOutCentCL1qual2 );
505 fOutputList->Add( fHOutMultV0MvsCL1qual2);
506 fOutputList->Add( fHOutMultV0MvsTRKqual2);
507 fOutputList->Add( fHOutMultTRKvsCL1qual2);
508 fOutputList->Add( fHOutQuality );
509 fOutputList->Add( fHOutVertex );
512 fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
514 PostData(1, fOutputList);
516 if (fPass==0) AliFatal("Which pass are you analyzing? You should set it via taskCentrality->SetPass(N)");
519 //________________________________________________________________________
520 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
522 // Execute analysis for current event:
523 if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
525 Float_t zncEnergy = 0.; // ZNC Energy
526 Float_t zpcEnergy = 0.; // ZPC Energy
527 Float_t znaEnergy = 0.; // ZNA Energy
528 Float_t zpaEnergy = 0.; // ZPA Energy
529 Float_t zem1Energy = 0.; // ZEM1 Energy
530 Float_t zem2Energy = 0.; // ZEM2 Energy
531 Bool_t zdcEnergyCal = kFALSE; // if zdc is calibrated (in pass2)
533 Int_t nTracks = 0; // no. tracks
534 Int_t nTracklets = 0; // no. tracklets
535 Int_t nClusters[6] = {0}; // no. clusters on 6 ITS layers
536 Int_t nChips[2]; // no. chips on 2 SPD layers
537 Float_t spdCorr =0; // corrected spd2 multiplicity
539 Float_t multV0A = 0; // multiplicity from V0 reco side A
540 Float_t multV0C = 0; // multiplicity from V0 reco side C
541 Short_t multV0AOnline = 0; // multiplicity from V0 reco side A
542 Short_t multV0COnline = 0; // multiplicity from V0 reco side C
543 Float_t v0Corr = 0; // corrected V0 multiplicity (used for MC)
545 Float_t multFMDA = 0; // multiplicity from FMD on detector A
546 Float_t multFMDC = 0; // multiplicity from FMD on detector C
548 Float_t zvtx =0; // z-vertex SPD
549 Int_t zvtxNcont =0; // contributors to z-vertex SPD
552 AliCentrality *esdCent = 0;
554 if(fAnalysisInput.CompareTo("ESD")==0){
556 AliVEvent* event = InputEvent();
557 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
559 AliError("No ESD Event");
565 if (SetupRun(esd)<0) {
566 AliError("Centrality File not available for this run");
570 esdCent = esd->GetCentrality();
573 AliESDVZERO* esdV0 = esd->GetVZEROData();
574 multV0A=esdV0->GetMTotV0A();
575 multV0C=esdV0->GetMTotV0C();
576 v0Corr = multV0A+multV0C;
578 multV0AOnline=esdV0->GetTriggerChargeA();
579 multV0COnline=esdV0->GetTriggerChargeC();
583 const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
584 zvtx = vtxESD->GetZ();
585 zvtxNcont = vtxESD->GetNContributors();
587 // ***** CB info (tracklets, clusters, chips)
588 //nTracks = event->GetNumberOfTracks();
589 nTracks = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
591 const AliMultiplicity *mult = esd->GetMultiplicity();
593 nTracklets = mult->GetNumberOfTracklets();
595 for(Int_t ilay=0; ilay<6; ilay++){
596 nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
599 for(Int_t ilay=0; ilay<2; ilay++){
600 nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
603 spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);
606 AliESDFMD *fmd = esd->GetFMDData();
607 Float_t totalMultA = 0;
608 Float_t totalMultC = 0;
609 const Float_t fFMDLowCut = 0.4;
611 for(UShort_t det=1;det<=3;det++) {
612 Int_t nRings = (det==1 ? 1 : 2);
613 for (UShort_t ir = 0; ir < nRings; ir++) {
614 Char_t ring = (ir == 0 ? 'I' : 'O');
615 UShort_t nsec = (ir == 0 ? 20 : 40);
616 UShort_t nstr = (ir == 0 ? 512 : 256);
617 for(UShort_t sec =0; sec < nsec; sec++) {
618 for(UShort_t strip = 0; strip < nstr; strip++) {
620 Float_t fmdMult = fmd->Multiplicity(det,ring,sec,strip);
621 if(fmdMult == 0 || fmdMult == AliESDFMD::kInvalidMult) continue;
623 Float_t nParticles=0;
625 if(fmdMult > fFMDLowCut) {
629 if (det<3) totalMultA = totalMultA + nParticles;
630 else totalMultC = totalMultC + nParticles;
636 multFMDA = totalMultA;
637 multFMDC = totalMultC;
640 AliESDZDC *esdZDC = esd->GetESDZDC();
641 zdcEnergyCal = esdZDC->AliESDZDC::TestBit(AliESDZDC::kEnergyCalibratedSignal);
643 zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
644 zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
645 znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
646 zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
648 zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
649 zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
650 znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
651 zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
653 zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
654 zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
658 else if(fAnalysisInput.CompareTo("AOD")==0){
659 //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
661 printf(" AOD analysis not yet implemented!!!\n\n");
665 // ***** Scaling for MC
668 v0Corr = Short_t((multV0A+multV0C) * fV0MScaleFactorMC);
670 // ***** Scaling for Data
672 v0Corr = Short_t(v0Corr / fV0MScaleFactor);
673 spdCorr = spdCorr / fSPDScaleFactor;
674 nTracks = Int_t(nTracks / fTPCScaleFactor);
677 // ***** Centrality Selection
678 if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
679 if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
680 if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
681 if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
682 if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
683 if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
685 if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
686 if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
687 if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin(zem1Energy+zem2Energy,zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
695 if (TMath::Abs(zvtx)>fZVCut || zvtxNcont<1) fQuality += 1;
697 // ***** outliers, skip in case of MC input
700 if (IsOutlierV0MSPD(spdCorr, v0Corr, int(fCentV0M))) fQuality += 2;
702 if (IsOutlierV0MTPC(nTracks, v0Corr, int(fCentV0M))) fQuality += 4;
704 if (IsOutlierV0MZDC((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
705 (zdcEnergyCal==kFALSE) ) fQuality += 8;
706 if (IsOutlierV0MZDCECal((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
707 (zdcEnergyCal==kTRUE) ) fQuality += 8;
715 esdCent->SetQuality(fQuality);
716 esdCent->SetCentralityV0M(fCentV0M);
717 esdCent->SetCentralityFMD(fCentFMD);
718 esdCent->SetCentralityTRK(fCentTRK);
719 esdCent->SetCentralityTKL(fCentTKL);
720 esdCent->SetCentralityCL0(fCentCL0);
721 esdCent->SetCentralityCL1(fCentCL1);
722 esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
723 esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
724 esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
727 fHOutQuality->Fill(fQuality);
728 fHOutVertex->Fill(zvtx);
731 fHOutCentV0M->Fill(fCentV0M);
732 fHOutCentFMD->Fill(fCentFMD);
733 fHOutCentTRK->Fill(fCentTRK);
734 fHOutCentTKL->Fill(fCentTKL);
735 fHOutCentCL0->Fill(fCentCL0);
736 fHOutCentCL1->Fill(fCentCL1);
737 fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
738 fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
739 fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
740 fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentCL1);
741 fHOutCentV0MvsCentTRK->Fill(fCentV0M,fCentTRK);
742 fHOutCentTRKvsCentCL1->Fill(fCentTRK,fCentCL1);
743 fHOutCentV0MvsCentZDC->Fill(fCentV0M,fCentZEMvsZDC);
744 fHOutMultV0M->Fill(multV0A+multV0C);
745 fHOutMultV0O->Fill(multV0AOnline+multV0COnline);
746 fHOutMultFMD->Fill(multFMDA+multFMDC);
747 fHOutMultTRK->Fill(nTracks);
748 fHOutMultTKL->Fill(nTracklets);
749 fHOutMultCL0->Fill(nClusters[0]);
750 fHOutMultCL1->Fill(spdCorr);
751 fHOutMultV0MvsZDN->Fill(v0Corr,(zncEnergy+znaEnergy));
752 fHOutMultZEMvsZDN->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy));
753 fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
754 fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
755 fHOutMultZEMvsZDCw->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy),fCentV0M);
756 fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
757 fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
758 fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
759 fHOutMultV0MvsV0O->Fill(v0Corr,(multV0AOnline+multV0COnline));
760 fHOutMultV0OvsCL1->Fill((multV0AOnline+multV0COnline),spdCorr);
761 fHOutMultV0OvsTRK->Fill((multV0AOnline+multV0COnline),nTracks);
762 } else if (fQuality%2 == 0) {
763 fHOutCentV0Mqual1->Fill(fCentV0M);
764 fHOutCentTRKqual1->Fill(fCentTRK);
765 fHOutCentCL1qual1->Fill(fCentCL1);
766 fHOutMultV0MvsCL1qual1->Fill(v0Corr,spdCorr);
767 fHOutMultV0MvsTRKqual1->Fill(v0Corr,nTracks);
768 fHOutMultTRKvsCL1qual1->Fill(nTracks,spdCorr);
770 fHOutCentV0Mqual2->Fill(fCentV0M);
771 fHOutCentTRKqual2->Fill(fCentTRK);
772 fHOutCentCL1qual2->Fill(fCentCL1);
773 fHOutMultV0MvsCL1qual2->Fill(v0Corr,spdCorr);
774 fHOutMultV0MvsTRKqual2->Fill(v0Corr,nTracks);
775 fHOutMultTRKvsCL1qual2->Fill(nTracks,spdCorr);
778 PostData(1, fOutputList);
780 //________________________________________________________________________
781 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
783 // Terminate analysis
785 //________________________________________________________________________
786 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* const esd)
788 // Setup files for run
793 // check if something to be done
794 if (fCurrentRun == esd->GetRunNumber())
797 fCurrentRun = esd->GetRunNumber();
799 TString fileName =(Form("%s/COMMON/CENTRALITY/data/centrality.root", AliAnalysisManager::GetOADBPath()));
800 AliInfo(Form("Setup Centrality Selection for run %d with file %s\n",fCurrentRun,fileName.Data()));
802 AliOADBContainer *con = new AliOADBContainer("OADB");
803 con->InitFromFile(fileName,"Centrality");
805 AliOADBCentrality* centOADB = 0;
806 centOADB = (AliOADBCentrality*)(con->GetObject(fCurrentRun));
808 AliWarning(Form("Centrality OADB does not exist for run %d, using Default \n",fCurrentRun ));
809 centOADB = (AliOADBCentrality*)(con->GetDefaultObject("oadbDefault"));
813 fUseScaling = centOADB->UseScaling();
814 fUseCleaning = centOADB->UseCleaning();
817 fZVCut = centOADB->ZVCut();
818 fOutliersCut = centOADB->OutliersCut();
821 fHtempV0M = centOADB->V0hist();
822 fHtempTRK = centOADB->TPChist();
823 fHtempCL1 = centOADB->SPDhist();
824 fHtempZEMvsZDC = centOADB->ZEMvsZDChist();
826 TString path = gSystem->ExpandPathName(fileName.Data());
827 if (!fHtempV0M) AliWarning(Form("Calibration for V0M does not exist in %s", path.Data()));
828 if (!fHtempTRK) AliWarning(Form("Calibration for TRK does not exist in %s", path.Data()));
829 if (!fHtempCL1) AliWarning(Form("Calibration for CL1 does not exist in %s", path.Data()));
830 if (!fHtempZEMvsZDC) AliWarning(Form("Calibration for ZEMvsZDC does not exist in %s", path.Data()));
833 fV0MScaleFactor = centOADB->V0MScaleFactor();
834 fSPDScaleFactor = centOADB->SPDScaleFactor();
835 fTPCScaleFactor = centOADB->TPCScaleFactor();
836 fV0MScaleFactorMC = centOADB->V0MScaleFactor();
838 // outliers parameters
839 fV0MSPDOutlierPar0 = centOADB->V0MSPDOutlierPar0();
840 fV0MSPDOutlierPar1 = centOADB->V0MSPDOutlierPar1();
841 fV0MTPCOutlierPar0 = centOADB->V0MTPCOutlierPar0();
842 fV0MTPCOutlierPar1 = centOADB->V0MTPCOutlierPar1();
844 fV0MSPDSigmaOutlierPar0 = centOADB->V0MSPDSigmaOutlierPar0();
845 fV0MSPDSigmaOutlierPar1 = centOADB->V0MSPDSigmaOutlierPar1();
846 fV0MSPDSigmaOutlierPar2 = centOADB->V0MSPDSigmaOutlierPar2();
847 fV0MTPCSigmaOutlierPar0 = centOADB->V0MTPCSigmaOutlierPar0();
848 fV0MTPCSigmaOutlierPar1 = centOADB->V0MTPCSigmaOutlierPar1();
849 fV0MTPCSigmaOutlierPar2 = centOADB->V0MTPCSigmaOutlierPar2();
851 fV0MZDCOutlierPar0 = centOADB->V0MZDCOutlierPar0();
852 fV0MZDCOutlierPar1 = centOADB->V0MZDCOutlierPar1();
853 fV0MZDCEcalOutlierPar0 = centOADB->V0MZDCEcalOutlierPar0();
854 fV0MZDCEcalOutlierPar1 = centOADB->V0MZDCEcalOutlierPar1();
863 //________________________________________________________________________
864 Bool_t AliCentralitySelectionTask::IsOutlierV0MSPD(Float_t spd, Float_t v0, Int_t cent) const
867 Float_t val = fV0MSPDOutlierPar0 + fV0MSPDOutlierPar1 * v0;
868 Float_t spdSigma = fV0MSPDSigmaOutlierPar0 + fV0MSPDSigmaOutlierPar1*cent + fV0MSPDSigmaOutlierPar2*cent*cent;
869 if ( TMath::Abs(spd-val) > fOutliersCut*spdSigma )
875 //________________________________________________________________________
876 Bool_t AliCentralitySelectionTask::IsOutlierV0MTPC(Int_t tracks, Float_t v0, Int_t cent) const
879 Float_t val = fV0MTPCOutlierPar0 + fV0MTPCOutlierPar1 * v0;
880 Float_t tpcSigma = fV0MTPCSigmaOutlierPar0 + fV0MTPCSigmaOutlierPar1*cent + fV0MTPCSigmaOutlierPar2*cent*cent;
881 if ( TMath::Abs(tracks-val) > fOutliersCut*tpcSigma )
887 //________________________________________________________________________
888 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDC(Float_t zdc, Float_t v0) const
891 Float_t val = fV0MZDCOutlierPar0 + fV0MZDCOutlierPar1 * v0;
898 //________________________________________________________________________
899 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDCECal(Float_t zdc, Float_t v0) const
902 Float_t val = fV0MZDCEcalOutlierPar0 + fV0MZDCEcalOutlierPar1 * v0;