Fix in the last caall to CleanOwnPrimaryVertex
[u/mrichter/AliRoot.git] / ANALYSIS / AliCentralitySelectionTask.cxx
CommitLineData
a540a9d3 1/**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//*****************************************************
17// Class AliCentralitySelectionTask
18// Class to analyze determine centrality
19// author: Alberica Toia
20//*****************************************************
21
22#include "AliCentralitySelectionTask.h"
23
24#include <TTree.h>
25#include <TList.h>
26#include <TH1F.h>
27#include <TH2F.h>
28#include <TF1.h>
29#include <TProfile.h>
30#include <TFile.h>
31#include <TObjString.h>
32#include <TString.h>
33#include <TCanvas.h>
34#include <TROOT.h>
35#include <TDirectory.h>
36#include <TSystem.h>
37#include <iostream>
38
39#include "AliAnalysisManager.h"
136f0a9b 40#include "AliHeader.h"
a540a9d3 41#include "AliVEvent.h"
42#include "AliESD.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"
a540a9d3 49#include "AliESDtrackCuts.h"
136f0a9b 50#include "AliESDVertex.h"
51#include "AliCentrality.h"
a540a9d3 52#include "AliMultiplicity.h"
53#include "AliAODHandler.h"
136f0a9b 54#include "AliAODHeader.h"
a540a9d3 55#include "AliAODEvent.h"
a540a9d3 56#include "AliAODVertex.h"
136f0a9b 57#include "AliAODVZERO.h"
58#include "AliAODTracklets.h"
a540a9d3 59#include "AliAODMCHeader.h"
a540a9d3 60#include "AliMCEventHandler.h"
136f0a9b 61#include "AliMCEvent.h"
62#include "AliAODMCParticle.h"
a540a9d3 63#include "AliMCParticle.h"
64#include "AliStack.h"
a540a9d3 65#include "AliAnalysisTaskSE.h"
66#include "AliGenEventHeader.h"
67#include "AliGenHijingEventHeader.h"
68#include "AliPhysicsSelectionTask.h"
69#include "AliPhysicsSelection.h"
70#include "AliBackgroundSelection.h"
71#include "AliESDUtils.h"
72
73ClassImp(AliCentralitySelectionTask)
74
75
76//________________________________________________________________________
77AliCentralitySelectionTask::AliCentralitySelectionTask():
78AliAnalysisTaskSE(),
a540a9d3 79 fAnalysisInput("ESD"),
80 fIsMCInput(kFALSE),
a121a512 81 fPass(0),
a540a9d3 82 fFile(0),
83 fFile2(0),
84 fCurrentRun(-1),
85 fRunNo(-1),
86 fLowRunN(0),
87 fHighRunN(0),
88 fUseScaling(0),
d015c169 89 fUseCleaning(0),
a540a9d3 90 fTrackCuts(0),
91 fZVCut(10),
92 fOutliersCut(5),
407bd85b 93 fQuality(999),
a540a9d3 94 fCentV0M(0),
95 fCentFMD(0),
96 fCentTRK(0),
97 fCentTKL(0),
98 fCentCL0(0),
99 fCentCL1(0),
100 fCentV0MvsFMD(0),
101 fCentTKLvsV0M(0),
102 fCentZEMvsZDC(0),
103 fHtempV0M(0),
104 fHtempFMD(0),
105 fHtempTRK(0),
106 fHtempTKL(0),
107 fHtempCL0(0),
108 fHtempCL1(0),
109 fHtempV0MvsFMD(0),
110 fHtempTKLvsV0M(0),
111 fHtempZEMvsZDC(0),
112 fOutputList(0),
113 fHOutCentV0M (0),
114 fHOutCentFMD (0),
115 fHOutCentTRK (0),
116 fHOutCentTKL (0),
117 fHOutCentCL0 (0),
118 fHOutCentCL1 (0),
119 fHOutCentV0MvsFMD(0),
120 fHOutCentTKLvsV0M(0),
121 fHOutCentZEMvsZDC(0),
9138a60c 122 fHOutCentV0MvsCentCL1(0),
123 fHOutCentV0MvsCentTRK(0),
124 fHOutCentTRKvsCentCL1(0),
980d6ad3 125 fHOutCentV0MvsCentZDC(0),
a540a9d3 126 fHOutMultV0M(0),
127 fHOutMultV0R(0),
128 fHOutMultFMD(0),
129 fHOutMultTRK(0),
130 fHOutMultTKL(0),
131 fHOutMultCL0(0),
132 fHOutMultCL1(0),
31200bdf 133 fHOutMultV0MvsZDN(0),
134 fHOutMultZEMvsZDN(0),
a540a9d3 135 fHOutMultV0MvsZDC(0),
136 fHOutMultZEMvsZDC(0),
39bbe018 137 fHOutMultZEMvsZDCw(0),
a540a9d3 138 fHOutMultV0MvsCL1(0),
139 fHOutMultV0MvsTRK(0),
140 fHOutMultTRKvsCL1(0),
04341115 141 fHOutCentV0Mqual1(0),
142 fHOutCentTRKqual1(0),
143 fHOutCentCL1qual1(0),
79d4544b 144 fHOutMultV0MvsCL1qual1(0),
145 fHOutMultV0MvsTRKqual1(0),
146 fHOutMultTRKvsCL1qual1(0),
04341115 147 fHOutCentV0Mqual2(0),
148 fHOutCentTRKqual2(0),
149 fHOutCentCL1qual2(0),
79d4544b 150 fHOutMultV0MvsCL1qual2(0),
151 fHOutMultV0MvsTRKqual2(0),
152 fHOutMultTRKvsCL1qual2(0),
a540a9d3 153 fHOutQuality(0),
154 fHOutVertex(0)
155{
156 // Default constructor
157 AliInfo("Centrality Selection enabled.");
609f601e 158 fLowRunN =136851;
159 fHighRunN=139517;
160
ff325947 161 for (Int_t i=0; i < 2667; i++) {
04341115 162 fV0MScaleFactor[i]=0.0;
163 fSPDScaleFactor[i]=0.0;
164 fTPCScaleFactor[i]=0.0;
165 fV0MScaleFactorMC[i]=0.0;
a540a9d3 166 }
167 fUseScaling=kTRUE;
d015c169 168 fUseCleaning=kTRUE;
cb15f505 169 fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO."
170 ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
a540a9d3 171}
172
173//________________________________________________________________________
174AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
175 AliAnalysisTaskSE(name),
a540a9d3 176 fAnalysisInput("ESD"),
177 fIsMCInput(kFALSE),
a121a512 178 fPass(0),
a540a9d3 179 fFile(0),
180 fFile2(0),
181 fCurrentRun(-1),
182 fRunNo(-1),
183 fLowRunN(0),
184 fHighRunN(0),
185 fUseScaling(0),
d015c169 186 fUseCleaning(0),
a540a9d3 187 fTrackCuts(0),
188 fZVCut(10),
189 fOutliersCut(5),
407bd85b 190 fQuality(999),
a540a9d3 191 fCentV0M(0),
192 fCentFMD(0),
193 fCentTRK(0),
194 fCentTKL(0),
195 fCentCL0(0),
196 fCentCL1(0),
197 fCentV0MvsFMD(0),
198 fCentTKLvsV0M(0),
199 fCentZEMvsZDC(0),
200 fHtempV0M(0),
201 fHtempFMD(0),
202 fHtempTRK(0),
203 fHtempTKL(0),
204 fHtempCL0(0),
205 fHtempCL1(0),
206 fHtempV0MvsFMD(0),
207 fHtempTKLvsV0M(0),
208 fHtempZEMvsZDC(0),
209 fOutputList(0),
210 fHOutCentV0M (0),
211 fHOutCentFMD (0),
212 fHOutCentTRK (0),
213 fHOutCentTKL (0),
214 fHOutCentCL0 (0),
215 fHOutCentCL1 (0),
216 fHOutCentV0MvsFMD(0),
217 fHOutCentTKLvsV0M(0),
218 fHOutCentZEMvsZDC(0),
9138a60c 219 fHOutCentV0MvsCentCL1(0),
220 fHOutCentV0MvsCentTRK(0),
221 fHOutCentTRKvsCentCL1(0),
980d6ad3 222 fHOutCentV0MvsCentZDC(0),
a540a9d3 223 fHOutMultV0M(0),
224 fHOutMultV0R(0),
225 fHOutMultFMD(0),
226 fHOutMultTRK(0),
227 fHOutMultTKL(0),
228 fHOutMultCL0(0),
229 fHOutMultCL1(0),
31200bdf 230 fHOutMultV0MvsZDN(0),
231 fHOutMultZEMvsZDN(0),
a540a9d3 232 fHOutMultV0MvsZDC(0),
233 fHOutMultZEMvsZDC(0),
39bbe018 234 fHOutMultZEMvsZDCw(0),
a540a9d3 235 fHOutMultV0MvsCL1(0),
236 fHOutMultV0MvsTRK(0),
237 fHOutMultTRKvsCL1(0),
04341115 238 fHOutCentV0Mqual1(0),
239 fHOutCentTRKqual1(0),
240 fHOutCentCL1qual1(0),
79d4544b 241 fHOutMultV0MvsCL1qual1(0),
242 fHOutMultV0MvsTRKqual1(0),
243 fHOutMultTRKvsCL1qual1(0),
04341115 244 fHOutCentV0Mqual2(0),
245 fHOutCentTRKqual2(0),
246 fHOutCentCL1qual2(0),
79d4544b 247 fHOutMultV0MvsCL1qual2(0),
248 fHOutMultV0MvsTRKqual2(0),
249 fHOutMultTRKvsCL1qual2(0),
a540a9d3 250 fHOutQuality(0),
251 fHOutVertex(0)
252{
609f601e 253 // Default constructor
254 fLowRunN =136851;
255 fHighRunN=139517;
256
a540a9d3 257 AliInfo("Centrality Selection enabled.");
258 DefineOutput(1, TList::Class());
ff325947 259 for (Int_t i=0; i<2667; i++) {
04341115 260 fV0MScaleFactor[i]=0.0;
261 fSPDScaleFactor[i]=0.0;
262 fTPCScaleFactor[i]=0.0;
263 fV0MScaleFactorMC[i]=0.0;
a540a9d3 264 }
265 fUseScaling=kTRUE;
d015c169 266 fUseCleaning=kTRUE;
cb15f505 267 fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO."
268 ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
a540a9d3 269}
270
271//________________________________________________________________________
272AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
273{
274 // Assignment operator
275 if (this!=&c) {
276 AliAnalysisTaskSE::operator=(c);
277 }
278 return *this;
279}
280
281//________________________________________________________________________
282AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
283 AliAnalysisTaskSE(ana),
a540a9d3 284 fAnalysisInput(ana.fDebug),
285 fIsMCInput(ana.fIsMCInput),
a121a512 286 fPass(ana.fPass),
a540a9d3 287 fFile(ana.fFile),
288 fFile2(ana.fFile2),
289 fCurrentRun(ana.fCurrentRun),
290 fRunNo(ana.fRunNo),
291 fLowRunN(ana.fLowRunN),
292 fHighRunN(ana.fHighRunN),
293 fUseScaling(ana.fUseScaling),
d015c169 294 fUseCleaning(ana.fUseCleaning),
a540a9d3 295 fTrackCuts(ana.fTrackCuts),
296 fZVCut(ana.fZVCut),
297 fOutliersCut(ana.fOutliersCut),
298 fQuality(ana.fQuality),
299 fCentV0M(ana.fCentV0M),
300 fCentFMD(ana.fCentFMD),
301 fCentTRK(ana.fCentTRK),
302 fCentTKL(ana.fCentTKL),
303 fCentCL0(ana.fCentCL0),
304 fCentCL1(ana.fCentCL1),
305 fCentV0MvsFMD(ana.fCentV0MvsFMD),
306 fCentTKLvsV0M(ana.fCentTKLvsV0M),
307 fCentZEMvsZDC(ana.fCentZEMvsZDC),
308 fHtempV0M(ana.fHtempV0M),
309 fHtempFMD(ana.fHtempFMD),
310 fHtempTRK(ana.fHtempTRK),
311 fHtempTKL(ana.fHtempTKL),
312 fHtempCL0(ana.fHtempCL0),
313 fHtempCL1(ana.fHtempCL1),
314 fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
315 fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
316 fHtempZEMvsZDC(ana.fHtempZEMvsZDC),
317 fOutputList(ana.fOutputList),
318 fHOutCentV0M (ana.fHOutCentV0M ),
319 fHOutCentFMD (ana.fHOutCentFMD ),
320 fHOutCentTRK (ana.fHOutCentTRK ),
321 fHOutCentTKL (ana.fHOutCentTKL ),
322 fHOutCentCL0 (ana.fHOutCentCL0 ),
323 fHOutCentCL1 (ana.fHOutCentCL1 ),
324 fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD),
325 fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M),
326 fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC),
9138a60c 327 fHOutCentV0MvsCentCL1(ana.fHOutCentV0MvsCentCL1),
328 fHOutCentV0MvsCentTRK(ana.fHOutCentV0MvsCentTRK),
329 fHOutCentTRKvsCentCL1(ana.fHOutCentTRKvsCentCL1),
980d6ad3 330 fHOutCentV0MvsCentZDC(ana.fHOutCentV0MvsCentZDC),
a540a9d3 331 fHOutMultV0M(ana.fHOutMultV0M),
332 fHOutMultV0R(ana.fHOutMultV0R),
333 fHOutMultFMD(ana.fHOutMultFMD),
334 fHOutMultTRK(ana.fHOutMultTRK),
335 fHOutMultTKL(ana.fHOutMultTKL),
336 fHOutMultCL0(ana.fHOutMultCL0),
337 fHOutMultCL1(ana.fHOutMultCL1),
31200bdf 338 fHOutMultV0MvsZDN(ana.fHOutMultV0MvsZDN),
339 fHOutMultZEMvsZDN(ana.fHOutMultZEMvsZDN),
a540a9d3 340 fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC),
341 fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),
39bbe018 342 fHOutMultZEMvsZDCw(ana.fHOutMultZEMvsZDCw),
a540a9d3 343 fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),
344 fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),
345 fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1),
04341115 346 fHOutCentV0Mqual1(ana.fHOutCentV0Mqual1),
347 fHOutCentTRKqual1(ana.fHOutCentTRKqual1),
348 fHOutCentCL1qual1(ana.fHOutCentCL1qual1),
79d4544b 349 fHOutMultV0MvsCL1qual1(ana.fHOutMultV0MvsCL1qual1),
350 fHOutMultV0MvsTRKqual1(ana.fHOutMultV0MvsTRKqual1),
351 fHOutMultTRKvsCL1qual1(ana.fHOutMultTRKvsCL1qual1),
04341115 352 fHOutCentV0Mqual2(ana.fHOutCentV0Mqual2),
353 fHOutCentTRKqual2(ana.fHOutCentTRKqual2),
354 fHOutCentCL1qual2(ana.fHOutCentCL1qual2),
79d4544b 355 fHOutMultV0MvsCL1qual2(ana.fHOutMultV0MvsCL1qual2),
356 fHOutMultV0MvsTRKqual2(ana.fHOutMultV0MvsTRKqual2),
357 fHOutMultTRKvsCL1qual2(ana.fHOutMultTRKvsCL1qual2),
a540a9d3 358 fHOutQuality(ana.fHOutQuality),
359 fHOutVertex(ana.fHOutVertex)
360{
361 // Copy Constructor
ff325947 362 for (Int_t i=0; i<2667; i++) {
04341115 363 fV0MScaleFactor[i]=0.0;
364 fSPDScaleFactor[i]=0.0;
365 fTPCScaleFactor[i]=0.0;
366 fV0MScaleFactorMC[i]=0.0;
7ee7a2df 367 }
368
a540a9d3 369}
370
371//________________________________________________________________________
372AliCentralitySelectionTask::~AliCentralitySelectionTask()
373{
374 // Destructor
375 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
376 if (fTrackCuts) delete fTrackCuts;
377}
378
379//________________________________________________________________________
380void AliCentralitySelectionTask::UserCreateOutputObjects()
381{
382 // Create the output containers
383 if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
384 AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
385
386 fOutputList = new TList();
387 fOutputList->SetOwner();
bab72d34 388 fHOutCentV0M = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",505,0,101);
389 fHOutCentFMD = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",505,0,101);
390 fHOutCentTRK = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",505,0,101);
391 fHOutCentTKL = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",505,0,101);
392 fHOutCentCL0 = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",505,0,101);
393 fHOutCentCL1 = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",505,0,101);
394 fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",505,0,101);
395 fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",505,0,101);
396 fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",505,0,101);
980d6ad3 397 fHOutCentV0MvsCentCL1= new TH2F("fHOutCentV0MvsCentCL1","fHOutCentV0MvsCentCL1; Cent V0; Cent SPD",505,0,101,505,0,101);
398 fHOutCentV0MvsCentTRK= new TH2F("fHOutCentV0MvsCentTRK","fHOutCentV0MvsCentTRK; Cent V0; Cent TPC",505,0,101,505,0,101);
399 fHOutCentTRKvsCentCL1= new TH2F("fHOutCentTRKvsCentCL1","fHOutCentTRKvsCentCL1; Cent TPC; Cent SPD",505,0,101,505,0,101);
400 fHOutCentV0MvsCentZDC= new TH2F("fHOutCentV0MvsCentZDC","fHOutCentV0MvsCentZDC; Cent V0; Cent ZDC",505,0,101,505,0,101);
a540a9d3 401
402 fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,25000);
2cf88888 403 fHOutMultV0R = new TH1F("fHOutMultV0R","fHOutMultV0R; Multiplicity V0",30000,0,30000);
a540a9d3 404 fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000);
405 fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000);
406 fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000);
407 fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000);
408 fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000);
31200bdf 409 fHOutMultV0MvsZDN = new TH2F("fHOutMultV0MvsZDN","fHOutMultV0MvsZDN; Multiplicity V0; Energy ZDC-N",500,0,25000,500,0,180000);
410 fHOutMultZEMvsZDN = new TH2F("fHOutMultZEMvsZDN","fHOutMultZEMvsZDN; Energy ZEM; Energy ZDC-N",500,0,2500,500,0,180000);
411 fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,25000,500,0,200000);
412 fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,200000);
39bbe018 413 fHOutMultZEMvsZDCw = new TH2F("fHOutMultZEMvsZDCw","fHOutMultZEMvsZDCw; Energy ZEM; Energy ZDC (weigthed with V0 percentile)",500,0,2500,500,0,200000);
a540a9d3 414 fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
415 fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
416 fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
417
04341115 418 fHOutCentV0Mqual1 = new TH1F("fHOutCentV0M_qual1","fHOutCentV0M_qual1; Centrality V0",505,0,101);
419 fHOutCentTRKqual1 = new TH1F("fHOutCentTRK_qual1","fHOutCentTRK_qual1; Centrality TPC",505,0,101);
420 fHOutCentCL1qual1 = new TH1F("fHOutCentCL1_qual1","fHOutCentCL1_qual1; Centrality SPD outer",505,0,101);
79d4544b 421 fHOutMultV0MvsCL1qual1 = new TH2F("fHOutMultV0MvsCL1_qual1","fHOutMultV0MvsCL1_qual1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
422 fHOutMultV0MvsTRKqual1 = new TH2F("fHOutMultV0MvsTRK_qual1","fHOutMultV0MvsTRK_qual1; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
423 fHOutMultTRKvsCL1qual1 = new TH2F("fHOutMultTRKvsCL1_qual1","fHOutMultTRKvsCL1_qual1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
a540a9d3 424
04341115 425 fHOutCentV0Mqual2 = new TH1F("fHOutCentV0M_qual2","fHOutCentV0M_qual2; Centrality V0",505,0,101);
426 fHOutCentTRKqual2 = new TH1F("fHOutCentTRK_qual2","fHOutCentTRK_qual2; Centrality TPC",505,0,101);
427 fHOutCentCL1qual2 = new TH1F("fHOutCentCL1_qual2","fHOutCentCL1_qual2; Centrality SPD outer",505,0,101);
79d4544b 428 fHOutMultV0MvsCL1qual2 = new TH2F("fHOutMultV0MvsCL1_qual2","fHOutMultV0MvsCL1_qual2; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
429 fHOutMultV0MvsTRKqual2 = new TH2F("fHOutMultV0MvsTRK_qual2","fHOutMultV0MvsTRK_qual2; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
430 fHOutMultTRKvsCL1qual2 = new TH2F("fHOutMultTRKvsCL1_qual2","fHOutMultTRKvsCL1_qual2; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
a540a9d3 431
407bd85b 432 fHOutQuality = new TH1F("fHOutQuality", "fHOutQuality", 100,-0.5,99.5);
a540a9d3 433 fHOutVertex = new TH1F("fHOutVertex", "fHOutVertex", 100,-20,20);
434
435 fOutputList->Add( fHOutCentV0M );
436 fOutputList->Add( fHOutCentFMD );
437 fOutputList->Add( fHOutCentTRK );
438 fOutputList->Add( fHOutCentTKL );
439 fOutputList->Add( fHOutCentCL0 );
440 fOutputList->Add( fHOutCentCL1 );
441 fOutputList->Add( fHOutCentV0MvsFMD);
442 fOutputList->Add( fHOutCentTKLvsV0M);
443 fOutputList->Add( fHOutCentZEMvsZDC);
9138a60c 444 fOutputList->Add( fHOutCentV0MvsCentCL1);
445 fOutputList->Add( fHOutCentV0MvsCentTRK);
446 fOutputList->Add( fHOutCentTRKvsCentCL1);
980d6ad3 447 fOutputList->Add( fHOutCentV0MvsCentZDC);
a540a9d3 448 fOutputList->Add( fHOutMultV0M);
449 fOutputList->Add( fHOutMultV0R);
450 fOutputList->Add( fHOutMultFMD);
451 fOutputList->Add( fHOutMultTRK);
452 fOutputList->Add( fHOutMultTKL);
453 fOutputList->Add( fHOutMultCL0);
454 fOutputList->Add( fHOutMultCL1);
31200bdf 455 fOutputList->Add( fHOutMultV0MvsZDN);
456 fOutputList->Add( fHOutMultZEMvsZDN);
a540a9d3 457 fOutputList->Add( fHOutMultV0MvsZDC);
458 fOutputList->Add( fHOutMultZEMvsZDC);
39bbe018 459 fOutputList->Add( fHOutMultZEMvsZDCw);
a540a9d3 460 fOutputList->Add( fHOutMultV0MvsCL1);
461 fOutputList->Add( fHOutMultV0MvsTRK);
462 fOutputList->Add( fHOutMultTRKvsCL1);
04341115 463 fOutputList->Add( fHOutCentV0Mqual1 );
464 fOutputList->Add( fHOutCentTRKqual1 );
465 fOutputList->Add( fHOutCentCL1qual1 );
79d4544b 466 fOutputList->Add( fHOutMultV0MvsCL1qual1);
467 fOutputList->Add( fHOutMultV0MvsTRKqual1);
468 fOutputList->Add( fHOutMultTRKvsCL1qual1);
04341115 469 fOutputList->Add( fHOutCentV0Mqual2 );
470 fOutputList->Add( fHOutCentTRKqual2 );
471 fOutputList->Add( fHOutCentCL1qual2 );
79d4544b 472 fOutputList->Add( fHOutMultV0MvsCL1qual2);
473 fOutputList->Add( fHOutMultV0MvsTRKqual2);
474 fOutputList->Add( fHOutMultTRKvsCL1qual2);
a540a9d3 475 fOutputList->Add( fHOutQuality );
476 fOutputList->Add( fHOutVertex );
477
478
479 fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
480
481 PostData(1, fOutputList);
482
483 MyInitScaleFactor();
484 if (fIsMCInput) MyInitScaleFactorMC();
afb1125b 485
486 if (fPass==0) AliFatal("Which pass are you analyzing? You should set it via taskCentrality->SetPass(N)");
a540a9d3 487}
488
489//________________________________________________________________________
490void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
491{
492 // Execute analysis for current event:
493 if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
494
495 Float_t zncEnergy = 0.; // ZNC Energy
496 Float_t zpcEnergy = 0.; // ZPC Energy
497 Float_t znaEnergy = 0.; // ZNA Energy
498 Float_t zpaEnergy = 0.; // ZPA Energy
499 Float_t zem1Energy = 0.; // ZEM1 Energy
500 Float_t zem2Energy = 0.; // ZEM2 Energy
cfead5a9 501 Bool_t zdcEnergyCal = kFALSE; // if zdc is calibrated (in pass2)
502
a540a9d3 503 Int_t nTracks = 0; // no. tracks
504 Int_t nTracklets = 0; // no. tracklets
505 Int_t nClusters[6] = {0}; // no. clusters on 6 ITS layers
506 Int_t nChips[2]; // no. chips on 2 SPD layers
507 Float_t spdCorr =0; // corrected spd2 multiplicity
508
509 Float_t multV0A = 0; // multiplicity from V0 reco side A
510 Float_t multV0C = 0; // multiplicity from V0 reco side C
511 Float_t multFMDA = 0; // multiplicity from FMD on detector A
512 Float_t multFMDC = 0; // multiplicity from FMD on detector C
513
514 Short_t v0Corr = 0; // corrected V0 multiplicity
515 Short_t v0CorrResc = 0; // corrected and rescaled V0 multiplicity
516
517 Float_t zvtx =0; // z-vertex SPD
407bd85b 518 Int_t zvtxNcont =0; // contributors to z-vertex SPD
136f0a9b 519
a540a9d3 520
521 AliCentrality *esdCent = 0;
522
523 if(fAnalysisInput.CompareTo("ESD")==0){
524
525 AliVEvent* event = InputEvent();
526 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
527 if (!esd) {
528 AliError("No ESD Event");
529 return;
530 }
cb15f505 531
532 LoadBranches();
a540a9d3 533
534 if (fRunNo<=0) {
dc6e8ace 535 if (SetupRun(esd)<0) {
536 AliError("Centrality File not available for this run");
537 return;
538 }
a540a9d3 539 }
540
541 esdCent = esd->GetCentrality();
542
543 // ***** V0 info
544 AliESDVZERO* esdV0 = esd->GetVZEROData();
545 multV0A=esdV0->GetMTotV0A();
546 multV0C=esdV0->GetMTotV0C();
547
548 Float_t v0CorrR;
549 v0Corr = (Short_t)AliESDUtils::GetCorrV0(esd,v0CorrR);
550 v0CorrResc = (Short_t)v0CorrR;
551
552 // ***** Vertex Info
553 const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
554 zvtx = vtxESD->GetZ();
407bd85b 555 zvtxNcont = vtxESD->GetNContributors();
a540a9d3 556
557 // ***** CB info (tracklets, clusters, chips)
558 //nTracks = event->GetNumberOfTracks();
559 nTracks = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
560
561 const AliMultiplicity *mult = esd->GetMultiplicity();
562
563 nTracklets = mult->GetNumberOfTracklets();
564
565 for(Int_t ilay=0; ilay<6; ilay++){
566 nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
567 }
568
569 for(Int_t ilay=0; ilay<2; ilay++){
570 nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
571 }
572
573 spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);
574
575 // ***** FMD info
576 AliESDFMD *fmd = esd->GetFMDData();
577 Float_t totalMultA = 0;
578 Float_t totalMultC = 0;
579 const Float_t fFMDLowCut = 0.4;
580
581 for(UShort_t det=1;det<=3;det++) {
582 Int_t nRings = (det==1 ? 1 : 2);
583 for (UShort_t ir = 0; ir < nRings; ir++) {
584 Char_t ring = (ir == 0 ? 'I' : 'O');
585 UShort_t nsec = (ir == 0 ? 20 : 40);
586 UShort_t nstr = (ir == 0 ? 512 : 256);
587 for(UShort_t sec =0; sec < nsec; sec++) {
588 for(UShort_t strip = 0; strip < nstr; strip++) {
589
04341115 590 Float_t fmdMult = fmd->Multiplicity(det,ring,sec,strip);
591 if(fmdMult == 0 || fmdMult == AliESDFMD::kInvalidMult) continue;
a540a9d3 592
593 Float_t nParticles=0;
594
04341115 595 if(fmdMult > fFMDLowCut) {
a540a9d3 596 nParticles = 1.;
597 }
598
599 if (det<3) totalMultA = totalMultA + nParticles;
600 else totalMultC = totalMultC + nParticles;
601
602 }
603 }
604 }
605 }
606 multFMDA = totalMultA;
607 multFMDC = totalMultC;
608
609 // ***** ZDC info
610 AliESDZDC *esdZDC = esd->GetESDZDC();
cfead5a9 611 zdcEnergyCal = esdZDC->AliESDZDC::TestBit(AliESDZDC::kEnergyCalibratedSignal);
612 if (zdcEnergyCal) {
613 zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
614 zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
615 znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
616 zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
617 } else {
618 zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
619 zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
620 znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
621 zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
622 }
a540a9d3 623 zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
624 zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
cfead5a9 625
a540a9d3 626 }
136f0a9b 627
a540a9d3 628 else if(fAnalysisInput.CompareTo("AOD")==0){
629 //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
630 // to be implemented
631 printf(" AOD analysis not yet implemented!!!\n\n");
632 return;
633 }
634
a540a9d3 635 // ***** Scaling
a121a512 636 // ***** Scaling for pass2
637 if (fPass==2) {
638 fUseScaling=kFALSE;
a121a512 639 }
d015c169 640 // ***** Scaling for MC
641 if (fIsMCInput) {
642 fUseScaling=kFALSE;
de5b7aa5 643 fUseCleaning=kFALSE;
04341115 644 Float_t tempScalefactorV0M = MyGetScaleFactorMC(fCurrentRun);
645 v0Corr = Short_t((multV0A+multV0C) * tempScalefactorV0M);
d015c169 646 }
647 // ***** Scaling for Data
a540a9d3 648 if (fUseScaling) {
04341115 649 Float_t tempScalefactorV0M = MyGetScaleFactor(fCurrentRun,0);
650 Float_t tempScalefactorSPD = MyGetScaleFactor(fCurrentRun,1);
651 Float_t tempScalefactorTPC = MyGetScaleFactor(fCurrentRun,2);
652 v0Corr = Short_t(v0Corr / tempScalefactorV0M);
653 spdCorr = spdCorr / tempScalefactorSPD;
654 nTracks = Int_t(nTracks / tempScalefactorTPC);
a540a9d3 655 }
a540a9d3 656
657 // ***** Centrality Selection
658 if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
659 if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
660 if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
661 if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
662 if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
663 if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
664
665 if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
666 if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
7ee7a2df 667 if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin(zem1Energy+zem2Energy,zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
a540a9d3 668
669 // ***** Cleaning
d015c169 670 if (fUseCleaning) {
671 fQuality=0;
672 fZVCut=10;
673 fOutliersCut=6;
674
675 // ***** vertex
407bd85b 676 if (TMath::Abs(zvtx)>fZVCut || zvtxNcont<1) fQuality += 1;
d015c169 677
678 // ***** outliers
679 // **** V0 vs SPD
680 if (IsOutlierV0MSPD(spdCorr, v0Corr, int(fCentV0M))) fQuality += 2;
681 // ***** V0 vs TPC
682 if (IsOutlierV0MTPC(nTracks, v0Corr, int(fCentV0M))) fQuality += 4;
683 // ***** V0 vs ZDC
684 if (IsOutlierV0MZDC((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
685 (zdcEnergyCal==kFALSE) && !(fIsMCInput)) fQuality += 8;
686 if (IsOutlierV0MZDCECal((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
687 ((zdcEnergyCal==kTRUE) || (fIsMCInput))) fQuality += 8;
688 } else {
689 fQuality = 0;
690 }
a540a9d3 691
d015c169 692
a540a9d3 693 if (esdCent) {
694 esdCent->SetQuality(fQuality);
695 esdCent->SetCentralityV0M(fCentV0M);
696 esdCent->SetCentralityFMD(fCentFMD);
697 esdCent->SetCentralityTRK(fCentTRK);
698 esdCent->SetCentralityTKL(fCentTKL);
699 esdCent->SetCentralityCL0(fCentCL0);
700 esdCent->SetCentralityCL1(fCentCL1);
701 esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
702 esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
703 esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
704 }
705
706 fHOutQuality->Fill(fQuality);
707 fHOutVertex->Fill(zvtx);
a540a9d3 708
709 if (fQuality==0) {
710 fHOutCentV0M->Fill(fCentV0M);
711 fHOutCentFMD->Fill(fCentFMD);
712 fHOutCentTRK->Fill(fCentTRK);
713 fHOutCentTKL->Fill(fCentTKL);
714 fHOutCentCL0->Fill(fCentCL0);
715 fHOutCentCL1->Fill(fCentCL1);
716 fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
717 fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
718 fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
9138a60c 719 fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentCL1);
720 fHOutCentV0MvsCentTRK->Fill(fCentV0M,fCentTRK);
721 fHOutCentTRKvsCentCL1->Fill(fCentTRK,fCentCL1);
39bbe018 722 fHOutCentV0MvsCentZDC->Fill(fCentV0M,fCentZEMvsZDC);
a540a9d3 723 fHOutMultV0M->Fill(v0Corr);
724 fHOutMultV0R->Fill(multV0A+multV0C);
725 fHOutMultFMD->Fill((multFMDA+multFMDC));
726 fHOutMultTRK->Fill(nTracks);
727 fHOutMultTKL->Fill(nTracklets);
728 fHOutMultCL0->Fill(nClusters[0]);
729 fHOutMultCL1->Fill(spdCorr);
79d4544b 730 fHOutMultV0MvsZDN->Fill(v0Corr,(zncEnergy+znaEnergy));
731 fHOutMultZEMvsZDN->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy));
732 fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
733 fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
39bbe018 734 fHOutMultZEMvsZDCw->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy),fCentV0M);
79d4544b 735 fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
736 fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
a540a9d3 737 fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
79d4544b 738 } else if (fQuality%2 == 0) {
04341115 739 fHOutCentV0Mqual1->Fill(fCentV0M);
740 fHOutCentTRKqual1->Fill(fCentTRK);
741 fHOutCentCL1qual1->Fill(fCentCL1);
79d4544b 742 fHOutMultV0MvsCL1qual1->Fill(v0Corr,spdCorr);
743 fHOutMultV0MvsTRKqual1->Fill(v0Corr,nTracks);
744 fHOutMultTRKvsCL1qual1->Fill(nTracks,spdCorr);
a540a9d3 745 } else {
04341115 746 fHOutCentV0Mqual2->Fill(fCentV0M);
747 fHOutCentTRKqual2->Fill(fCentTRK);
748 fHOutCentCL1qual2->Fill(fCentCL1);
79d4544b 749 fHOutMultV0MvsCL1qual2->Fill(v0Corr,spdCorr);
750 fHOutMultV0MvsTRKqual2->Fill(v0Corr,nTracks);
751 fHOutMultTRKvsCL1qual2->Fill(nTracks,spdCorr);
a540a9d3 752 }
753
754 PostData(1, fOutputList);
755}
756
757//________________________________________________________________________
758void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename)
759{
760 // Read centrality histograms
761 TDirectory *owd = gDirectory;
762 // Check if the file is present
763 TString path = gSystem->ExpandPathName(fCentfilename.Data());
764 if (gSystem->AccessPathName(path)) {
765 AliError(Form("File %s does not exist", path.Data()));
766 return;
767 }
768 fFile = TFile::Open(fCentfilename);
769 owd->cd();
770 fHtempV0M = (TH1F*) (fFile->Get("hmultV0_percentile"));
771 fHtempFMD = (TH1F*) (fFile->Get("hmultFMD_percentile"));
772 fHtempTRK = (TH1F*) (fFile->Get("hNtracks_percentile"));
773 fHtempTKL = (TH1F*) (fFile->Get("hNtracklets_percentile"));
774 fHtempCL0 = (TH1F*) (fFile->Get("hNclusters0_percentile"));
775 fHtempCL1 = (TH1F*) (fFile->Get("hNclusters1_percentile"));
776
777 if (!fHtempV0M) AliWarning(Form("Calibration for V0M does not exist in %s", path.Data()));
778 if (!fHtempFMD) AliWarning(Form("Calibration for FMD does not exist in %s", path.Data()));
779 if (!fHtempTRK) AliWarning(Form("Calibration for TRK does not exist in %s", path.Data()));
780 if (!fHtempTKL) AliWarning(Form("Calibration for TKL does not exist in %s", path.Data()));
781 if (!fHtempCL0) AliWarning(Form("Calibration for CL0 does not exist in %s", path.Data()));
782 if (!fHtempCL1) AliWarning(Form("Calibration for CL1 does not exist in %s", path.Data()));
783
784 owd->cd();
785}
786
787//________________________________________________________________________
788void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2)
789{
790 // Read centrality histograms
791 TDirectory *owd = gDirectory;
792 TString path = gSystem->ExpandPathName(fCentfilename2.Data());
793 if (gSystem->AccessPathName(path)) {
794 AliError(Form("File %s does not exist", path.Data()));
795 return;
796 }
797 fFile2 = TFile::Open(fCentfilename2);
798 owd->cd();
799 fHtempV0MvsFMD = (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
800 fHtempTKLvsV0M = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
7ee7a2df 801 fHtempZEMvsZDC = (TH2F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
a540a9d3 802
803 if (!fHtempV0MvsFMD) AliWarning(Form("Calibration for V0MvsFMD does not exist in %s", path.Data()));
804 if (!fHtempTKLvsV0M) AliWarning(Form("Calibration for TKLvsV0M does not exist in %s", path.Data()));
805 if (!fHtempZEMvsZDC) AliWarning(Form("Calibration for ZEMvsZDC does not exist in %s", path.Data()));
806
807 owd->cd();
808}
809
810//________________________________________________________________________
811void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
812{
813 // Terminate analysis
814 if (fFile && fFile->IsOpen())
815 fFile->Close();
816 if (fFile2 && fFile2->IsOpen())
817 fFile2->Close();
818}
819//________________________________________________________________________
f690bf48 820Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* const esd)
a540a9d3 821{
822 // Setup files for run
823
824 if (!esd)
825 return -1;
826
827 // check if something to be done
828 if (fCurrentRun == esd->GetRunNumber())
829 return 0;
830 else
831 fCurrentRun = esd->GetRunNumber();
832
833 AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
834
835 // CHANGE HERE FOR RUN RANGES
db1b4c23 836 switch (fPass) {
837 case 1:
838 if ( fCurrentRun <= 137165 ) fRunNo = 137161;
839 else fRunNo = 137366;
840 break;
841 case 2:
407bd85b 842 if ( fCurrentRun >= 136851 && fCurrentRun <= 137165 ) fRunNo = 137161;
0b622d40 843 else if ( fCurrentRun >= 137230 && fCurrentRun <= 137531 ) fRunNo = 137366;
844 else if ( fCurrentRun >= 137539 && fCurrentRun <= 137848 ) fRunNo = 137722;
980d6ad3 845 else if ( fCurrentRun >= 138125 && fCurrentRun <= 138154 ) fRunNo = 138150;
79d4544b 846 else if ( fCurrentRun >= 138190 && fCurrentRun <= 138275 ) fRunNo = 138200;
980d6ad3 847 else fRunNo = 139172;
db1b4c23 848 break;
849 default:
850 AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
851 return -1;
852 }
853 TString fileName =(Form("%s/COMMON/CENTRALITY/data/pass%d/AliCentralityBy1D_%d.root", AliAnalysisManager::GetOADBPath(), fPass, fRunNo));
854 TString fileName2=(Form("%s/COMMON/CENTRALITY/data/pass%d/AliCentralityByFunction_%d.root", AliAnalysisManager::GetOADBPath(), fPass, fRunNo));
a540a9d3 855
a121a512 856 AliInfo(Form("Centrality Selection for run %d and pass %d is initialized with %s", fCurrentRun, fPass, fileName.Data()));
a540a9d3 857 ReadCentralityHistos(fileName.Data());
858 ReadCentralityHistos2(fileName2.Data());
859 if (!fFile && !fFile2) {
74b089bb 860 AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
a540a9d3 861 return -1;
862 }
863 return 0;
864}
865
866//________________________________________________________________________
f690bf48 867Bool_t AliCentralitySelectionTask::IsOutlierV0MSPD(Float_t spd, Float_t v0, Int_t cent) const
a540a9d3 868{
f690bf48 869// Clean outliers
136f0a9b 870 Float_t val= -0.579579 + 0.273949 * v0;
871 Float_t spdSigma[101]={216.473, 194.377, 190.432, 185.019, 179.787, 175.99, 170.695, 167.276, 162.229, 159.016, 155.592, 151.909, 148.967, 146.367, 142.752, 139.718, 136.838, 134.133, 131.374, 128.616, 126.195, 123.49, 120.76, 118.368, 115.671, 113.354, 110.986, 108.761, 105.816, 103.373, 101.525, 99.3941, 96.8891, 95.1021, 92.728, 90.6434, 88.6393, 86.4624, 84.19, 82.3937, 80.4803, 78.5633, 76.5256, 74.7059, 73.0146, 70.9321, 69.1671, 67.2923, 65.5859, 63.9027, 61.7826, 60.1855, 58.6887, 56.8523, 55.1695, 53.422, 51.9776, 50.3506, 48.5304, 47.1517, 45.4786, 43.9914, 42.4579, 41.0889, 39.6201, 38.3046, 36.8624, 35.3697, 33.9223, 32.4637, 30.859, 29.9707, 28.764, 27.5677, 26.3437, 24.8656, 23.5253, 22.9878, 21.5172, 20.4645, 19.7614, 19.0363, 18.0875, 17.3675, 16.7313, 16.4101, 15.8235, 15.4795, 14.897, 14.578, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 22.5965, 22.5965};
a540a9d3 872
04341115 873 if ( TMath::Abs(spd-val) > fOutliersCut*spdSigma[cent] )
a540a9d3 874 return kTRUE;
875 else
876 return kFALSE;
877}
878
879//________________________________________________________________________
f690bf48 880Bool_t AliCentralitySelectionTask::IsOutlierV0MTPC(Int_t tracks, Float_t v0, Int_t cent) const
a540a9d3 881{
f690bf48 882// Clean outliers
136f0a9b 883 Float_t val = -1.03873 + 0.125691 * v0;
884 Float_t tpcSigma[101]={105.1, 95.7594, 94.5777, 91.7369, 89.655, 87.9469, 85.2747, 83.8137, 81.5663, 79.9871, 78.3491, 76.664, 75.1352, 73.9838, 72.0122, 70.7627, 69.1832, 67.8707, 66.3974, 65.0222, 64.1065, 62.7248, 61.498, 60.11, 58.7405, 57.681, 56.2979, 55.3153, 53.983, 52.6813, 51.5321, 50.5531, 49.3309, 48.2269, 47.1665, 46.1924, 45.1143, 44.0394, 42.8571, 41.8471, 40.9594, 39.8298, 38.8128, 37.8845, 36.8788, 35.9896, 34.9477, 34.1293, 33.1267, 32.161, 31.1132, 30.3521, 29.5237, 28.5992, 27.6283, 26.8885, 26.0719, 25.1334, 24.2125, 23.5376, 22.5383, 21.8396, 21.0384, 20.2942, 19.4746, 18.5545, 18.26, 17.2997, 16.7109, 15.9324, 15.2755, 14.764, 14.2358, 13.6728, 13.2606, 12.846, 12.2636, 11.722, 11.1491, 10.5952, 10.235, 9.87852, 9.48122, 9.17571, 8.90995, 8.81675, 8.62595, 8.59371, 8.52936, 8.58058, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 14.2584, 14.2584};
a540a9d3 885
04341115 886 if ( TMath::Abs(tracks-val) > fOutliersCut*tpcSigma[cent] )
a540a9d3 887 return kTRUE;
888 else
889 return kFALSE;
890}
891
892//________________________________________________________________________
6abef102 893Bool_t AliCentralitySelectionTask::IsOutlierV0MZDC(Float_t zdc, Float_t v0) const
a540a9d3 894{
f690bf48 895// Clean outliers
6abef102 896 Float_t val = 235000. - 9.5 * v0;
897 if (zdc > val)
898 return kTRUE;
899 else
a540a9d3 900 return kFALSE;
901}
902
9138a60c 903//________________________________________________________________________
f690bf48 904Bool_t AliCentralitySelectionTask::IsOutlierV0MZDCECal(Float_t /*zdc*/, Float_t /*v0*/) const
9138a60c 905{
f690bf48 906// Clean outliers
9138a60c 907 return kFALSE;
908}
909
a540a9d3 910//________________________________________________________________________
f690bf48 911Float_t AliCentralitySelectionTask::MyGetScaleFactor(Int_t runnumber, Int_t flag) const
a540a9d3 912{
f690bf48 913// Get scaling factor
a540a9d3 914 if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
915 cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
916 return 0.0;
917 }
918
919 Float_t scalefactor=0.0;
920 if (flag==0)
04341115 921 scalefactor = fV0MScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
a540a9d3 922 else if (flag==1)
04341115 923 scalefactor = fSPDScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
a540a9d3 924 else if (flag==2)
04341115 925 scalefactor = fTPCScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
a540a9d3 926
927 return scalefactor;
928
929}
930
931//________________________________________________________________________
f690bf48 932Float_t AliCentralitySelectionTask::MyGetScaleFactorMC(Int_t runnumber) const
a540a9d3 933{
f690bf48 934// Get MC scaling factor
a540a9d3 935 if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
936 cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
937 return 0.0;
938 }
939
04341115 940 Float_t scalefactor= fV0MScaleFactorMC[runnumber - fLowRunN]; // subtracting reference offset index
a540a9d3 941 return scalefactor;
942
943}
944
945//________________________________________________________________________
946void AliCentralitySelectionTask::MyInitScaleFactor ()
947{
f690bf48 948// Initialize the scaling factors
04341115 949 for (int i=0; i<=(fHighRunN-fLowRunN); i++) fV0MScaleFactor[i] = 0.0;
950 for (int i=0; i<=(fHighRunN-fLowRunN); i++) fSPDScaleFactor[i] = 0.0;
951 for (int i=0; i<=(fHighRunN-fLowRunN); i++) fTPCScaleFactor[i] = 0.0;
a540a9d3 952
953 // scale factors determined from <V0 charge> on a run-by-run basis
04341115 954 fV0MScaleFactor[310] = 1.;
955 fV0MScaleFactor[311] = 1.;
956 fV0MScaleFactor[514] = 1.0046;
957 fV0MScaleFactor[515] = 0.983535;
958 fV0MScaleFactor[579] = 0.988185;
959 fV0MScaleFactor[580] = 0.983351;
960 fV0MScaleFactor[581] = 0.989013;
961 fV0MScaleFactor[583] = 0.990056;
962 fV0MScaleFactor[588] = 0.974438;
963 fV0MScaleFactor[589] = 0.981572;
964 fV0MScaleFactor[590] = 0.989316;
965 fV0MScaleFactor[592] = 0.98345;
966 fV0MScaleFactor[688] = 0.993647;
967 fV0MScaleFactor[690] = 0.994758;
968 fV0MScaleFactor[693] = 0.989569;
969 fV0MScaleFactor[698] = 0.993119;
970 fV0MScaleFactor[744] = 0.989583;
971 fV0MScaleFactor[757] = 0.990377;
972 fV0MScaleFactor[787] = 0.990176;
973 fV0MScaleFactor[788] = 0.98723;
974 fV0MScaleFactor[834] = 1.00403;
975 fV0MScaleFactor[835] = 0.994376;
976 fV0MScaleFactor[840] = 0.99759;
977 fV0MScaleFactor[842] = 1.01031;
978 fV0MScaleFactor[900] = 0.996216;
979 fV0MScaleFactor[901] = 0.994205;
980 fV0MScaleFactor[992] = 0.998479;
981 fV0MScaleFactor[997] = 1.00078;
982 fV0MScaleFactor[1299] = 1.00515;
983 fV0MScaleFactor[1303] = 1.00094;
984 fV0MScaleFactor[1339] = 0.986596;
985 fV0MScaleFactor[1346] = 0.972226;
986 fV0MScaleFactor[1349] = 0.960358;
987 fV0MScaleFactor[1350] = 0.970023;
988 fV0MScaleFactor[1374] = 1.00575;
989 fV0MScaleFactor[1545] = 1.00471;
990 fV0MScaleFactor[1587] = 1.00611;
991 fV0MScaleFactor[1588] = 1.00976;
992 fV0MScaleFactor[1618] = 1.00771;
993 fV0MScaleFactor[1682] = 1.01622;
994 fV0MScaleFactor[1727] = 1.01305;
995 fV0MScaleFactor[1728] = 1.00388;
996 fV0MScaleFactor[1731] = 1.00673;
997 fV0MScaleFactor[1732] = 1.00916;
998 fV0MScaleFactor[1770] = 1.0092;
999 fV0MScaleFactor[1773] = 1.00728;
1000 fV0MScaleFactor[1786] = 1.01655;
1001 fV0MScaleFactor[1787] = 1.00672;
1002 fV0MScaleFactor[1801] = 0.983339;
1003 fV0MScaleFactor[1802] = 1.00754;
1004 fV0MScaleFactor[1811] = 1.00608;
1005 fV0MScaleFactor[1815] = 1.01227;
1006 fV0MScaleFactor[1879] = 0.99907;
1007 fV0MScaleFactor[1881] = 0.995696;
1008 fV0MScaleFactor[2186] = 1.00559;
1009 fV0MScaleFactor[2187] = 1.00631;
1010 fV0MScaleFactor[2254] = 1.01512;
1011 fV0MScaleFactor[2256] = 0.998727;
1012 fV0MScaleFactor[2321] = 1.00701;
1013 fSPDScaleFactor[310] = 1.00211;
1014 fSPDScaleFactor[311] = 1.00067;
1015 fSPDScaleFactor[514] = 1.02268;
1016 fSPDScaleFactor[515] = 0.994902;
1017 fSPDScaleFactor[579] = 1.00215;
1018 fSPDScaleFactor[580] = 0.993421;
1019 fSPDScaleFactor[581] = 1.00129;
1020 fSPDScaleFactor[583] = 1.00242;
1021 fSPDScaleFactor[588] = 0.984762;
1022 fSPDScaleFactor[589] = 0.994355;
1023 fSPDScaleFactor[590] = 1.00073;
1024 fSPDScaleFactor[592] = 0.995889;
1025 fSPDScaleFactor[688] = 0.994532;
1026 fSPDScaleFactor[690] = 0.998307;
1027 fSPDScaleFactor[693] = 0.994052;
1028 fSPDScaleFactor[698] = 0.993224;
1029 fSPDScaleFactor[744] = 0.993279;
1030 fSPDScaleFactor[757] = 0.992494;
1031 fSPDScaleFactor[787] = 0.992678;
1032 fSPDScaleFactor[788] = 0.996563;
1033 fSPDScaleFactor[834] = 1.01116;
1034 fSPDScaleFactor[835] = 0.993108;
1035 fSPDScaleFactor[840] = 0.997574;
1036 fSPDScaleFactor[842] = 1.01829;
1037 fSPDScaleFactor[900] = 0.999438;
1038 fSPDScaleFactor[901] = 0.995849;
1039 fSPDScaleFactor[992] = 0.999227;
1040 fSPDScaleFactor[997] = 1.00575;
1041 fSPDScaleFactor[1299] = 0.99877;
1042 fSPDScaleFactor[1303] = 0.999682;
1043 fSPDScaleFactor[1339] = 0.978198;
1044 fSPDScaleFactor[1346] = 0.964178;
1045 fSPDScaleFactor[1349] = 0.959439;
1046 fSPDScaleFactor[1350] = 0.956945;
1047 fSPDScaleFactor[1374] = 0.994434;
1048 fSPDScaleFactor[1545] = 1.0016;
1049 fSPDScaleFactor[1587] = 1.00153;
1050 fSPDScaleFactor[1588] = 1.00698;
1051 fSPDScaleFactor[1618] = 1.00554;
1052 fSPDScaleFactor[1682] = 1.0123;
1053 fSPDScaleFactor[1727] = 1.011;
1054 fSPDScaleFactor[1728] = 1.00143;
1055 fSPDScaleFactor[1731] = 1.00486;
1056 fSPDScaleFactor[1732] = 1.00646;
1057 fSPDScaleFactor[1770] = 1.00515;
1058 fSPDScaleFactor[1773] = 1.00485;
1059 fSPDScaleFactor[1786] = 1.01215;
1060 fSPDScaleFactor[1787] = 1.00419;
1061 fSPDScaleFactor[1801] = 0.983327;
1062 fSPDScaleFactor[1802] = 1.00529;
1063 fSPDScaleFactor[1811] = 1.00367;
1064 fSPDScaleFactor[1815] = 1.01045;
1065 fSPDScaleFactor[1879] = 0.996374;
1066 fSPDScaleFactor[1881] = 0.988827;
1067 fSPDScaleFactor[2186] = 1.00354;
1068 fSPDScaleFactor[2187] = 1.00397;
1069 fSPDScaleFactor[2254] = 1.01138;
1070 fSPDScaleFactor[2256] = 0.996641;
1071 fSPDScaleFactor[2321] = 1.00357;
1072 fTPCScaleFactor[310] = 1.00434;
1073 fTPCScaleFactor[311] = 1.0056;
1074 fTPCScaleFactor[514] = 1.02185;
1075 fTPCScaleFactor[515] = 1.0036;
1076 fTPCScaleFactor[579] = 1.00607;
1077 fTPCScaleFactor[580] = 1.00183;
1078 fTPCScaleFactor[581] = 1.00693;
1079 fTPCScaleFactor[583] = 1.00746;
1080 fTPCScaleFactor[588] = 0.990524;
1081 fTPCScaleFactor[589] = 0.998582;
1082 fTPCScaleFactor[590] = 1.00618;
1083 fTPCScaleFactor[592] = 1.00088;
1084 fTPCScaleFactor[688] = 1.00598;
1085 fTPCScaleFactor[690] = 1.00658;
1086 fTPCScaleFactor[693] = 1.00144;
1087 fTPCScaleFactor[698] = 1.00279;
1088 fTPCScaleFactor[744] = 1.00122;
1089 fTPCScaleFactor[757] = 1.002;
1090 fTPCScaleFactor[787] = 0.997818;
1091 fTPCScaleFactor[788] = 0.994583;
1092 fTPCScaleFactor[834] = 1.01508;
1093 fTPCScaleFactor[835] = 1.00218;
1094 fTPCScaleFactor[840] = 1.00569;
1095 fTPCScaleFactor[842] = 1.01789;
1096 fTPCScaleFactor[900] = 1.00739;
1097 fTPCScaleFactor[901] = 1.00462;
1098 fTPCScaleFactor[992] = 1.00502;
1099 fTPCScaleFactor[997] = 1.00943;
1100 fTPCScaleFactor[1299] = 1.00438;
1101 fTPCScaleFactor[1303] = 0.996701;
1102 fTPCScaleFactor[1339] = 0.978641;
1103 fTPCScaleFactor[1346] = 0.968906;
1104 fTPCScaleFactor[1349] = 0.954311;
1105 fTPCScaleFactor[1350] = 0.958764;
1106 fTPCScaleFactor[1374] = 0.997899;
1107 fTPCScaleFactor[1545] = 0.992;
1108 fTPCScaleFactor[1587] = 0.992635;
1109 fTPCScaleFactor[1588] = 1.00207;
1110 fTPCScaleFactor[1618] = 1.00126;
1111 fTPCScaleFactor[1682] = 1.00324;
1112 fTPCScaleFactor[1727] = 1.00042;
1113 fTPCScaleFactor[1728] = 0.978881;
1114 fTPCScaleFactor[1731] = 0.999818;
1115 fTPCScaleFactor[1732] = 1.00109;
1116 fTPCScaleFactor[1770] = 0.99935;
1117 fTPCScaleFactor[1773] = 0.998531;
1118 fTPCScaleFactor[1786] = 0.999125;
1119 fTPCScaleFactor[1787] = 0.998479;
1120 fTPCScaleFactor[1801] = 0.9775;
1121 fTPCScaleFactor[1802] = 0.999095;
1122 fTPCScaleFactor[1811] = 0.998197;
1123 fTPCScaleFactor[1815] = 1.00413;
1124 fTPCScaleFactor[1879] = 0.990916;
1125 fTPCScaleFactor[1881] = 0.987241;
1126 fTPCScaleFactor[2186] = 1.00048;
1127 fTPCScaleFactor[2187] = 1.00057;
1128 fTPCScaleFactor[2254] = 1.00588;
1129 fTPCScaleFactor[2256] = 0.991503;
1130 fTPCScaleFactor[2321] = 1.00057;
609f601e 1131
a540a9d3 1132
1133 // set all missing values to the value of the run before it ....
bab72d34 1134 for (int i=0; i<=(fHighRunN-fLowRunN); i++) {
04341115 1135 if (fV0MScaleFactor[i] == 0.0) {
a540a9d3 1136 if (i==0) {
04341115 1137 fV0MScaleFactor[i] = 1.00;
a540a9d3 1138 } else {
1139 // search for last run number with non-zero value
1140 for (int j=i-1;j>=0;j--) {
04341115 1141 if (fV0MScaleFactor[j] != 0.0) {
1142 fV0MScaleFactor[i] = fV0MScaleFactor[j];
a540a9d3 1143 break;
1144 }
1145 }
1146 }
1147 }
1148 } // end loop over checking all run-numbers
1149
bab72d34 1150 for (int i=0; i<=(fHighRunN-fLowRunN); i++) {
04341115 1151 if (fSPDScaleFactor[i] == 0.0) {
a540a9d3 1152 if (i==0) {
04341115 1153 fSPDScaleFactor[i] = 1.00;
a540a9d3 1154 } else {
1155 for (int j=i-1;j>=0;j--) {
04341115 1156 if (fSPDScaleFactor[j] != 0.0) {
1157 fSPDScaleFactor[i] = fSPDScaleFactor[j];
a540a9d3 1158 break;
1159 }
1160 }
1161 }
1162 }
1163 }
1164
bab72d34 1165 for (int i=0; i<=(fHighRunN-fLowRunN); i++) {
04341115 1166 if (fTPCScaleFactor[i] == 0.0) {
a540a9d3 1167 if (i==0) {
04341115 1168 fTPCScaleFactor[i] = 1.00;
a540a9d3 1169 } else {
1170 for (int j=i-1;j>=0;j--) {
04341115 1171 if (fTPCScaleFactor[j] != 0.0) {
1172 fTPCScaleFactor[i] = fTPCScaleFactor[j];
a540a9d3 1173 break;
1174 }
1175 }
1176 }
1177 }
1178 }
1179
1180
04341115 1181 // for (int i=0; i<=(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << fV0MScaleFactor[i] << " for Run " << i+fLowRunN << endl;
a540a9d3 1182
1183 return;
1184
1185}
1186
1187
1188//________________________________________________________________________
1189void AliCentralitySelectionTask::MyInitScaleFactorMC()
1190{
f690bf48 1191// Initialize the MC scaling factors
04341115 1192 for (int i=0; i<(fHighRunN-fLowRunN); i++) fV0MScaleFactorMC[i] = 0.0;
a540a9d3 1193 // scale factors determined from <V0 charge> on a run-by-run basis
591abcd7 1194 switch (fPass) {
1195 case 1:
1196 fV0MScaleFactorMC[0] = 0.75108;
1197 break;
1198 case 2:
1199 fV0MScaleFactorMC[0] = 0.8;
1200 break;
1201 default:
1202 AliError(Form("Unknown reconstruction pass (%d)! Setting MC scale in V0 to 1", fPass));
1203 fV0MScaleFactorMC[0] = 1.0;
1204 break;
1205 }
1206
a540a9d3 1207 // set all missing values to the value of the run before it ....
bab72d34 1208 for (int i=0; i<=(fHighRunN-fLowRunN); i++) {
04341115 1209 if (fV0MScaleFactorMC[i] == 0.0) {
a540a9d3 1210 if (i==0) {
04341115 1211 fV0MScaleFactorMC[i] = 1.00;
a540a9d3 1212 } else {
1213 // search for last run number with non-zero value
1214 for (int j=i-1;j>=0;j--) {
04341115 1215 if (fV0MScaleFactorMC[j] != 0.0) {
1216 fV0MScaleFactorMC[i] = fV0MScaleFactorMC[j];
a540a9d3 1217 break;
1218 }
1219 }
1220 }
1221 }
1222 } // end loop over checking all run-numbers
1223
1224
04341115 1225 // for (int i=0; i<=(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << fV0MScaleFactorMC[i] << " for Run " << i+fLowRunN << endl;
a540a9d3 1226
1227 return;
1228
1229}
1230