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