]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ANALYSIS/AliCentralitySelectionTask.cxx
Optimisation
[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"
40#include "AliVEvent.h"
41#include "AliESD.h"
42#include "AliESDEvent.h"
43#include "AliESDHeader.h"
44#include "AliESDInputHandler.h"
45#include "AliESDZDC.h"
46#include "AliESDFMD.h"
47#include "AliESDVZERO.h"
48#include "AliCentrality.h"
49#include "AliESDtrackCuts.h"
50#include "AliMultiplicity.h"
51#include "AliAODHandler.h"
52#include "AliAODEvent.h"
53#include "AliESDVertex.h"
54#include "AliAODVertex.h"
55#include "AliAODMCHeader.h"
56#include "AliMCEvent.h"
57#include "AliMCEventHandler.h"
58#include "AliMCParticle.h"
59#include "AliStack.h"
60#include "AliHeader.h"
61#include "AliAODMCParticle.h"
62#include "AliAnalysisTaskSE.h"
63#include "AliGenEventHeader.h"
64#include "AliGenHijingEventHeader.h"
65#include "AliPhysicsSelectionTask.h"
66#include "AliPhysicsSelection.h"
67#include "AliBackgroundSelection.h"
68#include "AliESDUtils.h"
69
70ClassImp(AliCentralitySelectionTask)
71
72
73//________________________________________________________________________
74AliCentralitySelectionTask::AliCentralitySelectionTask():
75AliAnalysisTaskSE(),
76 fDebug(0),
77 fAnalysisInput("ESD"),
78 fIsMCInput(kFALSE),
79 fFile(0),
80 fFile2(0),
81 fCurrentRun(-1),
82 fRunNo(-1),
83 fLowRunN(0),
84 fHighRunN(0),
85 fUseScaling(0),
d015c169 86 fUseCleaning(0),
a540a9d3 87 fTrackCuts(0),
88 fZVCut(10),
89 fOutliersCut(5),
90 fQuality(0),
91 fCentV0M(0),
92 fCentFMD(0),
93 fCentTRK(0),
94 fCentTKL(0),
95 fCentCL0(0),
96 fCentCL1(0),
97 fCentV0MvsFMD(0),
98 fCentTKLvsV0M(0),
99 fCentZEMvsZDC(0),
100 fHtempV0M(0),
101 fHtempFMD(0),
102 fHtempTRK(0),
103 fHtempTKL(0),
104 fHtempCL0(0),
105 fHtempCL1(0),
106 fHtempV0MvsFMD(0),
107 fHtempTKLvsV0M(0),
108 fHtempZEMvsZDC(0),
109 fOutputList(0),
110 fHOutCentV0M (0),
111 fHOutCentFMD (0),
112 fHOutCentTRK (0),
113 fHOutCentTKL (0),
114 fHOutCentCL0 (0),
115 fHOutCentCL1 (0),
116 fHOutCentV0MvsFMD(0),
117 fHOutCentTKLvsV0M(0),
118 fHOutCentZEMvsZDC(0),
9138a60c 119 fHOutCentV0MvsCentCL1(0),
120 fHOutCentV0MvsCentTRK(0),
121 fHOutCentTRKvsCentCL1(0),
a540a9d3 122 fHOutMultV0M(0),
123 fHOutMultV0R(0),
124 fHOutMultFMD(0),
125 fHOutMultTRK(0),
126 fHOutMultTKL(0),
127 fHOutMultCL0(0),
128 fHOutMultCL1(0),
31200bdf 129 fHOutMultV0MvsZDN(0),
130 fHOutMultZEMvsZDN(0),
a540a9d3 131 fHOutMultV0MvsZDC(0),
132 fHOutMultZEMvsZDC(0),
133 fHOutMultV0MvsCL1(0),
134 fHOutMultV0MvsTRK(0),
135 fHOutMultTRKvsCL1(0),
136 fHOutCentV0M_qual1(0),
137 fHOutCentTRK_qual1(0),
138 fHOutCentCL1_qual1(0),
139 fHOutCentV0M_qual2(0),
140 fHOutCentTRK_qual2(0),
141 fHOutCentCL1_qual2(0),
142 fHOutQuality(0),
143 fHOutVertex(0)
144{
145 // Default constructor
146 AliInfo("Centrality Selection enabled.");
609f601e 147 fLowRunN =136851;
148 fHighRunN=139517;
149
ff325947 150 for (Int_t i=0; i < 2667; i++) {
a540a9d3 151 V0MScaleFactor[i]=0.0;
152 SPDScaleFactor[i]=0.0;
153 TPCScaleFactor[i]=0.0;
154 V0MScaleFactorMC[i]=0.0;
155 }
156 fUseScaling=kTRUE;
d015c169 157 fUseCleaning=kTRUE;
cb15f505 158 fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO."
159 ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
a540a9d3 160}
161
162//________________________________________________________________________
163AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
164 AliAnalysisTaskSE(name),
165 fDebug(0),
166 fAnalysisInput("ESD"),
167 fIsMCInput(kFALSE),
168 fFile(0),
169 fFile2(0),
170 fCurrentRun(-1),
171 fRunNo(-1),
172 fLowRunN(0),
173 fHighRunN(0),
174 fUseScaling(0),
d015c169 175 fUseCleaning(0),
a540a9d3 176 fTrackCuts(0),
177 fZVCut(10),
178 fOutliersCut(5),
179 fQuality(0),
180 fCentV0M(0),
181 fCentFMD(0),
182 fCentTRK(0),
183 fCentTKL(0),
184 fCentCL0(0),
185 fCentCL1(0),
186 fCentV0MvsFMD(0),
187 fCentTKLvsV0M(0),
188 fCentZEMvsZDC(0),
189 fHtempV0M(0),
190 fHtempFMD(0),
191 fHtempTRK(0),
192 fHtempTKL(0),
193 fHtempCL0(0),
194 fHtempCL1(0),
195 fHtempV0MvsFMD(0),
196 fHtempTKLvsV0M(0),
197 fHtempZEMvsZDC(0),
198 fOutputList(0),
199 fHOutCentV0M (0),
200 fHOutCentFMD (0),
201 fHOutCentTRK (0),
202 fHOutCentTKL (0),
203 fHOutCentCL0 (0),
204 fHOutCentCL1 (0),
205 fHOutCentV0MvsFMD(0),
206 fHOutCentTKLvsV0M(0),
207 fHOutCentZEMvsZDC(0),
9138a60c 208 fHOutCentV0MvsCentCL1(0),
209 fHOutCentV0MvsCentTRK(0),
210 fHOutCentTRKvsCentCL1(0),
a540a9d3 211 fHOutMultV0M(0),
212 fHOutMultV0R(0),
213 fHOutMultFMD(0),
214 fHOutMultTRK(0),
215 fHOutMultTKL(0),
216 fHOutMultCL0(0),
217 fHOutMultCL1(0),
31200bdf 218 fHOutMultV0MvsZDN(0),
219 fHOutMultZEMvsZDN(0),
a540a9d3 220 fHOutMultV0MvsZDC(0),
221 fHOutMultZEMvsZDC(0),
222 fHOutMultV0MvsCL1(0),
223 fHOutMultV0MvsTRK(0),
224 fHOutMultTRKvsCL1(0),
225 fHOutCentV0M_qual1(0),
226 fHOutCentTRK_qual1(0),
227 fHOutCentCL1_qual1(0),
228 fHOutCentV0M_qual2(0),
229 fHOutCentTRK_qual2(0),
230 fHOutCentCL1_qual2(0),
231 fHOutQuality(0),
232 fHOutVertex(0)
233{
609f601e 234 // Default constructor
235 fLowRunN =136851;
236 fHighRunN=139517;
237
a540a9d3 238 AliInfo("Centrality Selection enabled.");
239 DefineOutput(1, TList::Class());
ff325947 240 for (Int_t i=0; i<2667; i++) {
a540a9d3 241 V0MScaleFactor[i]=0.0;
242 SPDScaleFactor[i]=0.0;
243 TPCScaleFactor[i]=0.0;
244 V0MScaleFactorMC[i]=0.0;
245 }
246 fUseScaling=kTRUE;
d015c169 247 fUseCleaning=kTRUE;
cb15f505 248 fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO."
249 ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
a540a9d3 250}
251
252//________________________________________________________________________
253AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
254{
255 // Assignment operator
256 if (this!=&c) {
257 AliAnalysisTaskSE::operator=(c);
258 }
259 return *this;
260}
261
262//________________________________________________________________________
263AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
264 AliAnalysisTaskSE(ana),
265 fDebug(ana.fDebug),
266 fAnalysisInput(ana.fDebug),
267 fIsMCInput(ana.fIsMCInput),
268 fFile(ana.fFile),
269 fFile2(ana.fFile2),
270 fCurrentRun(ana.fCurrentRun),
271 fRunNo(ana.fRunNo),
272 fLowRunN(ana.fLowRunN),
273 fHighRunN(ana.fHighRunN),
274 fUseScaling(ana.fUseScaling),
d015c169 275 fUseCleaning(ana.fUseCleaning),
a540a9d3 276 fTrackCuts(ana.fTrackCuts),
277 fZVCut(ana.fZVCut),
278 fOutliersCut(ana.fOutliersCut),
279 fQuality(ana.fQuality),
280 fCentV0M(ana.fCentV0M),
281 fCentFMD(ana.fCentFMD),
282 fCentTRK(ana.fCentTRK),
283 fCentTKL(ana.fCentTKL),
284 fCentCL0(ana.fCentCL0),
285 fCentCL1(ana.fCentCL1),
286 fCentV0MvsFMD(ana.fCentV0MvsFMD),
287 fCentTKLvsV0M(ana.fCentTKLvsV0M),
288 fCentZEMvsZDC(ana.fCentZEMvsZDC),
289 fHtempV0M(ana.fHtempV0M),
290 fHtempFMD(ana.fHtempFMD),
291 fHtempTRK(ana.fHtempTRK),
292 fHtempTKL(ana.fHtempTKL),
293 fHtempCL0(ana.fHtempCL0),
294 fHtempCL1(ana.fHtempCL1),
295 fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
296 fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
297 fHtempZEMvsZDC(ana.fHtempZEMvsZDC),
298 fOutputList(ana.fOutputList),
299 fHOutCentV0M (ana.fHOutCentV0M ),
300 fHOutCentFMD (ana.fHOutCentFMD ),
301 fHOutCentTRK (ana.fHOutCentTRK ),
302 fHOutCentTKL (ana.fHOutCentTKL ),
303 fHOutCentCL0 (ana.fHOutCentCL0 ),
304 fHOutCentCL1 (ana.fHOutCentCL1 ),
305 fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD),
306 fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M),
307 fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC),
9138a60c 308 fHOutCentV0MvsCentCL1(ana.fHOutCentV0MvsCentCL1),
309 fHOutCentV0MvsCentTRK(ana.fHOutCentV0MvsCentTRK),
310 fHOutCentTRKvsCentCL1(ana.fHOutCentTRKvsCentCL1),
a540a9d3 311 fHOutMultV0M(ana.fHOutMultV0M),
312 fHOutMultV0R(ana.fHOutMultV0R),
313 fHOutMultFMD(ana.fHOutMultFMD),
314 fHOutMultTRK(ana.fHOutMultTRK),
315 fHOutMultTKL(ana.fHOutMultTKL),
316 fHOutMultCL0(ana.fHOutMultCL0),
317 fHOutMultCL1(ana.fHOutMultCL1),
31200bdf 318 fHOutMultV0MvsZDN(ana.fHOutMultV0MvsZDN),
319 fHOutMultZEMvsZDN(ana.fHOutMultZEMvsZDN),
a540a9d3 320 fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC),
321 fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),
322 fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),
323 fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),
324 fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1),
325 fHOutCentV0M_qual1(ana.fHOutCentV0M_qual1),
326 fHOutCentTRK_qual1(ana.fHOutCentTRK_qual1),
327 fHOutCentCL1_qual1(ana.fHOutCentCL1_qual1),
328 fHOutCentV0M_qual2(ana.fHOutCentV0M_qual2),
329 fHOutCentTRK_qual2(ana.fHOutCentTRK_qual2),
330 fHOutCentCL1_qual2(ana.fHOutCentCL1_qual2),
331 fHOutQuality(ana.fHOutQuality),
332 fHOutVertex(ana.fHOutVertex)
333{
334 // Copy Constructor
ff325947 335 for (Int_t i=0; i<2667; i++) {
7ee7a2df 336 V0MScaleFactor[i]=0.0;
337 SPDScaleFactor[i]=0.0;
338 TPCScaleFactor[i]=0.0;
339 V0MScaleFactorMC[i]=0.0;
340 }
341
a540a9d3 342}
343
344//________________________________________________________________________
345AliCentralitySelectionTask::~AliCentralitySelectionTask()
346{
347 // Destructor
348 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
349 if (fTrackCuts) delete fTrackCuts;
350}
351
352//________________________________________________________________________
353void AliCentralitySelectionTask::UserCreateOutputObjects()
354{
355 // Create the output containers
356 if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
357 AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
358
359 fOutputList = new TList();
360 fOutputList->SetOwner();
361 fHOutCentV0M = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",501,0,101);
362 fHOutCentFMD = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",501,0,101);
363 fHOutCentTRK = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",501,0,101);
364 fHOutCentTKL = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",501,0,101);
365 fHOutCentCL0 = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",501,0,101);
366 fHOutCentCL1 = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",501,0,101);
367 fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",501,0,101);
368 fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",501,0,101);
369 fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",501,0,101);
9138a60c 370 fHOutCentV0MvsCentCL1= new TH2F("fHOutCentV0MvsCentCL1","fHOutCentV0MvsCentCL1; Cent V0 vs Cent SPD",501,0,101,501,0,101);
371 fHOutCentV0MvsCentTRK= new TH2F("fHOutCentV0MvsCentTRK","fHOutCentV0MvsCentTRK; Cent V0 vs Cent TPC",501,0,101,501,0,101);
372 fHOutCentTRKvsCentCL1= new TH2F("fHOutCentTRKvsCentCL1","fHOutCentTRKvsCentCL1; Cent TPC vs Cent SPD",501,0,101,501,0,101);
a540a9d3 373
374 fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,25000);
2cf88888 375 fHOutMultV0R = new TH1F("fHOutMultV0R","fHOutMultV0R; Multiplicity V0",30000,0,30000);
a540a9d3 376 fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000);
377 fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000);
378 fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000);
379 fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000);
380 fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000);
31200bdf 381 fHOutMultV0MvsZDN = new TH2F("fHOutMultV0MvsZDN","fHOutMultV0MvsZDN; Multiplicity V0; Energy ZDC-N",500,0,25000,500,0,180000);
382 fHOutMultZEMvsZDN = new TH2F("fHOutMultZEMvsZDN","fHOutMultZEMvsZDN; Energy ZEM; Energy ZDC-N",500,0,2500,500,0,180000);
383 fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,25000,500,0,200000);
384 fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,200000);
a540a9d3 385 fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
386 fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
387 fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
388
389 fHOutCentV0M_qual1 = new TH1F("fHOutCentV0M_qual1","fHOutCentV0M_qual1; Centrality V0",501,0,100);
390 fHOutCentTRK_qual1 = new TH1F("fHOutCentTRK_qual1","fHOutCentTRK_qual1; Centrality TPC",501,0,100);
391 fHOutCentCL1_qual1 = new TH1F("fHOutCentCL1_qual1","fHOutCentCL1_qual1; Centrality SPD outer",501,0,100);
392
393 fHOutCentV0M_qual2 = new TH1F("fHOutCentV0M_qual2","fHOutCentV0M_qual2; Centrality V0",101,-0.1,100.1);
394 fHOutCentTRK_qual2 = new TH1F("fHOutCentTRK_qual2","fHOutCentTRK_qual2; Centrality TPC",101,-0.1,100.1);
395 fHOutCentCL1_qual2 = new TH1F("fHOutCentCL1_qual2","fHOutCentCL1_qual2; Centrality SPD outer",101,-0.1,100.1);
396
397 fHOutQuality = new TH1F("fHOutQuality", "fHOutQuality", 10,-0.5,9.5);
398 fHOutVertex = new TH1F("fHOutVertex", "fHOutVertex", 100,-20,20);
399
400 fOutputList->Add( fHOutCentV0M );
401 fOutputList->Add( fHOutCentFMD );
402 fOutputList->Add( fHOutCentTRK );
403 fOutputList->Add( fHOutCentTKL );
404 fOutputList->Add( fHOutCentCL0 );
405 fOutputList->Add( fHOutCentCL1 );
406 fOutputList->Add( fHOutCentV0MvsFMD);
407 fOutputList->Add( fHOutCentTKLvsV0M);
408 fOutputList->Add( fHOutCentZEMvsZDC);
9138a60c 409 fOutputList->Add( fHOutCentV0MvsCentCL1);
410 fOutputList->Add( fHOutCentV0MvsCentTRK);
411 fOutputList->Add( fHOutCentTRKvsCentCL1);
a540a9d3 412 fOutputList->Add( fHOutMultV0M);
413 fOutputList->Add( fHOutMultV0R);
414 fOutputList->Add( fHOutMultFMD);
415 fOutputList->Add( fHOutMultTRK);
416 fOutputList->Add( fHOutMultTKL);
417 fOutputList->Add( fHOutMultCL0);
418 fOutputList->Add( fHOutMultCL1);
31200bdf 419 fOutputList->Add( fHOutMultV0MvsZDN);
420 fOutputList->Add( fHOutMultZEMvsZDN);
a540a9d3 421 fOutputList->Add( fHOutMultV0MvsZDC);
422 fOutputList->Add( fHOutMultZEMvsZDC);
423 fOutputList->Add( fHOutMultV0MvsCL1);
424 fOutputList->Add( fHOutMultV0MvsTRK);
425 fOutputList->Add( fHOutMultTRKvsCL1);
426 fOutputList->Add( fHOutCentV0M_qual1 );
427 fOutputList->Add( fHOutCentTRK_qual1 );
428 fOutputList->Add( fHOutCentCL1_qual1 );
429 fOutputList->Add( fHOutCentV0M_qual2 );
430 fOutputList->Add( fHOutCentTRK_qual2 );
431 fOutputList->Add( fHOutCentCL1_qual2 );
432 fOutputList->Add( fHOutQuality );
433 fOutputList->Add( fHOutVertex );
434
435
436 fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
437
438 PostData(1, fOutputList);
439
440 MyInitScaleFactor();
441 if (fIsMCInput) MyInitScaleFactorMC();
442}
443
444//________________________________________________________________________
445void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
446{
447 // Execute analysis for current event:
448 if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
449
450 Float_t zncEnergy = 0.; // ZNC Energy
451 Float_t zpcEnergy = 0.; // ZPC Energy
452 Float_t znaEnergy = 0.; // ZNA Energy
453 Float_t zpaEnergy = 0.; // ZPA Energy
454 Float_t zem1Energy = 0.; // ZEM1 Energy
455 Float_t zem2Energy = 0.; // ZEM2 Energy
cfead5a9 456 Bool_t zdcEnergyCal = kFALSE; // if zdc is calibrated (in pass2)
457
a540a9d3 458 Int_t nTracks = 0; // no. tracks
459 Int_t nTracklets = 0; // no. tracklets
460 Int_t nClusters[6] = {0}; // no. clusters on 6 ITS layers
461 Int_t nChips[2]; // no. chips on 2 SPD layers
462 Float_t spdCorr =0; // corrected spd2 multiplicity
463
464 Float_t multV0A = 0; // multiplicity from V0 reco side A
465 Float_t multV0C = 0; // multiplicity from V0 reco side C
466 Float_t multFMDA = 0; // multiplicity from FMD on detector A
467 Float_t multFMDC = 0; // multiplicity from FMD on detector C
468
469 Short_t v0Corr = 0; // corrected V0 multiplicity
470 Short_t v0CorrResc = 0; // corrected and rescaled V0 multiplicity
471
472 Float_t zvtx =0; // z-vertex SPD
473
474 AliCentrality *esdCent = 0;
475
476 if(fAnalysisInput.CompareTo("ESD")==0){
477
478 AliVEvent* event = InputEvent();
479 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
480 if (!esd) {
481 AliError("No ESD Event");
482 return;
483 }
cb15f505 484
485 LoadBranches();
a540a9d3 486
487 if (fRunNo<=0) {
488 if (SetupRun(esd)<0)
489 AliFatal("Centrality File not available for this run");
490 }
491
492 esdCent = esd->GetCentrality();
493
494 // ***** V0 info
495 AliESDVZERO* esdV0 = esd->GetVZEROData();
496 multV0A=esdV0->GetMTotV0A();
497 multV0C=esdV0->GetMTotV0C();
498
499 Float_t v0CorrR;
500 v0Corr = (Short_t)AliESDUtils::GetCorrV0(esd,v0CorrR);
501 v0CorrResc = (Short_t)v0CorrR;
502
503 // ***** Vertex Info
504 const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
505 zvtx = vtxESD->GetZ();
506
507 // ***** CB info (tracklets, clusters, chips)
508 //nTracks = event->GetNumberOfTracks();
509 nTracks = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
510
511 const AliMultiplicity *mult = esd->GetMultiplicity();
512
513 nTracklets = mult->GetNumberOfTracklets();
514
515 for(Int_t ilay=0; ilay<6; ilay++){
516 nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
517 }
518
519 for(Int_t ilay=0; ilay<2; ilay++){
520 nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
521 }
522
523 spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);
524
525 // ***** FMD info
526 AliESDFMD *fmd = esd->GetFMDData();
527 Float_t totalMultA = 0;
528 Float_t totalMultC = 0;
529 const Float_t fFMDLowCut = 0.4;
530
531 for(UShort_t det=1;det<=3;det++) {
532 Int_t nRings = (det==1 ? 1 : 2);
533 for (UShort_t ir = 0; ir < nRings; ir++) {
534 Char_t ring = (ir == 0 ? 'I' : 'O');
535 UShort_t nsec = (ir == 0 ? 20 : 40);
536 UShort_t nstr = (ir == 0 ? 512 : 256);
537 for(UShort_t sec =0; sec < nsec; sec++) {
538 for(UShort_t strip = 0; strip < nstr; strip++) {
539
540 Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
541 if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
542
543 Float_t nParticles=0;
544
545 if(FMDmult > fFMDLowCut) {
546 nParticles = 1.;
547 }
548
549 if (det<3) totalMultA = totalMultA + nParticles;
550 else totalMultC = totalMultC + nParticles;
551
552 }
553 }
554 }
555 }
556 multFMDA = totalMultA;
557 multFMDC = totalMultC;
558
559 // ***** ZDC info
560 AliESDZDC *esdZDC = esd->GetESDZDC();
cfead5a9 561 zdcEnergyCal = esdZDC->AliESDZDC::TestBit(AliESDZDC::kEnergyCalibratedSignal);
562 if (zdcEnergyCal) {
563 zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
564 zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
565 znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
566 zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
567 } else {
568 zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
569 zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
570 znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
571 zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
572 }
a540a9d3 573 zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
574 zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
cfead5a9 575
a540a9d3 576 }
577 else if(fAnalysisInput.CompareTo("AOD")==0){
578 //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
579 // to be implemented
580 printf(" AOD analysis not yet implemented!!!\n\n");
581 return;
582 }
583
584
585 // ***** Scaling
d015c169 586 // ***** Scaling for MC
587 if (fIsMCInput) {
588 fUseScaling=kFALSE;
589 Float_t temp_scalefactorV0M = MyGetScaleFactorMC(fCurrentRun);
590 v0Corr = Short_t((multV0A+multV0C) * temp_scalefactorV0M);
591 }
592 // ***** Scaling for Data
a540a9d3 593 if (fUseScaling) {
594 Float_t temp_scalefactorV0M = MyGetScaleFactor(fCurrentRun,0);
595 Float_t temp_scalefactorSPD = MyGetScaleFactor(fCurrentRun,1);
596 Float_t temp_scalefactorTPC = MyGetScaleFactor(fCurrentRun,2);
597 v0Corr = Short_t(v0Corr / temp_scalefactorV0M);
598 spdCorr = spdCorr / temp_scalefactorSPD;
599 nTracks = Int_t(nTracks / temp_scalefactorTPC);
600 }
a540a9d3 601
602 // ***** Centrality Selection
603 if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
604 if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
605 if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
606 if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
607 if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
608 if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
609
610 if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
611 if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
7ee7a2df 612 if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin(zem1Energy+zem2Energy,zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
a540a9d3 613
614 // ***** Cleaning
d015c169 615 if (fUseCleaning) {
616 fQuality=0;
617 fZVCut=10;
618 fOutliersCut=6;
619
620 // ***** vertex
621 if (TMath::Abs(zvtx)>fZVCut) fQuality += 1;
622
623 // ***** outliers
624 // **** V0 vs SPD
625 if (IsOutlierV0MSPD(spdCorr, v0Corr, int(fCentV0M))) fQuality += 2;
626 // ***** V0 vs TPC
627 if (IsOutlierV0MTPC(nTracks, v0Corr, int(fCentV0M))) fQuality += 4;
628 // ***** V0 vs ZDC
629 if (IsOutlierV0MZDC((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
630 (zdcEnergyCal==kFALSE) && !(fIsMCInput)) fQuality += 8;
631 if (IsOutlierV0MZDCECal((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
632 ((zdcEnergyCal==kTRUE) || (fIsMCInput))) fQuality += 8;
633 } else {
634 fQuality = 0;
635 }
a540a9d3 636
d015c169 637
a540a9d3 638 if (esdCent) {
639 esdCent->SetQuality(fQuality);
640 esdCent->SetCentralityV0M(fCentV0M);
641 esdCent->SetCentralityFMD(fCentFMD);
642 esdCent->SetCentralityTRK(fCentTRK);
643 esdCent->SetCentralityTKL(fCentTKL);
644 esdCent->SetCentralityCL0(fCentCL0);
645 esdCent->SetCentralityCL1(fCentCL1);
646 esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
647 esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
648 esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
649 }
650
651 fHOutQuality->Fill(fQuality);
652 fHOutVertex->Fill(zvtx);
653
31200bdf 654 fHOutMultV0MvsZDN->Fill(v0Corr,(zncEnergy+znaEnergy));
655 fHOutMultZEMvsZDN->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy));
a540a9d3 656 fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
657 fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
658 fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
659 fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
660
661 if (fQuality==0) {
662 fHOutCentV0M->Fill(fCentV0M);
663 fHOutCentFMD->Fill(fCentFMD);
664 fHOutCentTRK->Fill(fCentTRK);
665 fHOutCentTKL->Fill(fCentTKL);
666 fHOutCentCL0->Fill(fCentCL0);
667 fHOutCentCL1->Fill(fCentCL1);
668 fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
669 fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
670 fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
9138a60c 671 fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentCL1);
672 fHOutCentV0MvsCentTRK->Fill(fCentV0M,fCentTRK);
673 fHOutCentTRKvsCentCL1->Fill(fCentTRK,fCentCL1);
a540a9d3 674 fHOutMultV0M->Fill(v0Corr);
675 fHOutMultV0R->Fill(multV0A+multV0C);
676 fHOutMultFMD->Fill((multFMDA+multFMDC));
677 fHOutMultTRK->Fill(nTracks);
678 fHOutMultTKL->Fill(nTracklets);
679 fHOutMultCL0->Fill(nClusters[0]);
680 fHOutMultCL1->Fill(spdCorr);
681 fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
682 } else if (fQuality ==1) {
683 fHOutCentV0M_qual1->Fill(fCentV0M);
684 fHOutCentTRK_qual1->Fill(fCentTRK);
685 fHOutCentCL1_qual1->Fill(fCentCL1);
686 } else {
687 fHOutCentV0M_qual2->Fill(fCentV0M);
688 fHOutCentTRK_qual2->Fill(fCentTRK);
689 fHOutCentCL1_qual2->Fill(fCentCL1);
690 }
691
692 PostData(1, fOutputList);
693}
694
695//________________________________________________________________________
696void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename)
697{
698 // Read centrality histograms
699 TDirectory *owd = gDirectory;
700 // Check if the file is present
701 TString path = gSystem->ExpandPathName(fCentfilename.Data());
702 if (gSystem->AccessPathName(path)) {
703 AliError(Form("File %s does not exist", path.Data()));
704 return;
705 }
706 fFile = TFile::Open(fCentfilename);
707 owd->cd();
708 fHtempV0M = (TH1F*) (fFile->Get("hmultV0_percentile"));
709 fHtempFMD = (TH1F*) (fFile->Get("hmultFMD_percentile"));
710 fHtempTRK = (TH1F*) (fFile->Get("hNtracks_percentile"));
711 fHtempTKL = (TH1F*) (fFile->Get("hNtracklets_percentile"));
712 fHtempCL0 = (TH1F*) (fFile->Get("hNclusters0_percentile"));
713 fHtempCL1 = (TH1F*) (fFile->Get("hNclusters1_percentile"));
714
715 if (!fHtempV0M) AliWarning(Form("Calibration for V0M does not exist in %s", path.Data()));
716 if (!fHtempFMD) AliWarning(Form("Calibration for FMD does not exist in %s", path.Data()));
717 if (!fHtempTRK) AliWarning(Form("Calibration for TRK does not exist in %s", path.Data()));
718 if (!fHtempTKL) AliWarning(Form("Calibration for TKL does not exist in %s", path.Data()));
719 if (!fHtempCL0) AliWarning(Form("Calibration for CL0 does not exist in %s", path.Data()));
720 if (!fHtempCL1) AliWarning(Form("Calibration for CL1 does not exist in %s", path.Data()));
721
722 owd->cd();
723}
724
725//________________________________________________________________________
726void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2)
727{
728 // Read centrality histograms
729 TDirectory *owd = gDirectory;
730 TString path = gSystem->ExpandPathName(fCentfilename2.Data());
731 if (gSystem->AccessPathName(path)) {
732 AliError(Form("File %s does not exist", path.Data()));
733 return;
734 }
735 fFile2 = TFile::Open(fCentfilename2);
736 owd->cd();
737 fHtempV0MvsFMD = (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
738 fHtempTKLvsV0M = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
7ee7a2df 739 fHtempZEMvsZDC = (TH2F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
a540a9d3 740
741 if (!fHtempV0MvsFMD) AliWarning(Form("Calibration for V0MvsFMD does not exist in %s", path.Data()));
742 if (!fHtempTKLvsV0M) AliWarning(Form("Calibration for TKLvsV0M does not exist in %s", path.Data()));
743 if (!fHtempZEMvsZDC) AliWarning(Form("Calibration for ZEMvsZDC does not exist in %s", path.Data()));
744
745 owd->cd();
746}
747
748//________________________________________________________________________
749void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
750{
751 // Terminate analysis
752 if (fFile && fFile->IsOpen())
753 fFile->Close();
754 if (fFile2 && fFile2->IsOpen())
755 fFile2->Close();
756}
757//________________________________________________________________________
f690bf48 758Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* const esd)
a540a9d3 759{
760 // Setup files for run
761
762 if (!esd)
763 return -1;
764
765 // check if something to be done
766 if (fCurrentRun == esd->GetRunNumber())
767 return 0;
768 else
769 fCurrentRun = esd->GetRunNumber();
770
771 AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
772
773 // CHANGE HERE FOR RUN RANGES
774 if ( fCurrentRun <= 137165 ) fRunNo = 137161;
775 else fRunNo = 137366;
776 // CHANGE HERE FOR RUN RANGES
777
778 TString fileName(Form("%s/COMMON/CENTRALITY/data/AliCentralityBy1D_%d.root", AliAnalysisManager::GetOADBPath(), fRunNo));
779 TString fileName2(Form("%s/COMMON/CENTRALITY/data/AliCentralityByFunction_%d.root", AliAnalysisManager::GetOADBPath(), fRunNo));
780
781 AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));
782 ReadCentralityHistos(fileName.Data());
783 ReadCentralityHistos2(fileName2.Data());
784 if (!fFile && !fFile2) {
785 AliFatal(Form("Run %d not known to centrality selection!", fCurrentRun));
786 return -1;
787 }
788 return 0;
789}
790
791//________________________________________________________________________
f690bf48 792Bool_t AliCentralitySelectionTask::IsOutlierV0MSPD(Float_t spd, Float_t v0, Int_t cent) const
a540a9d3 793{
f690bf48 794// Clean outliers
f7292781 795 Float_t val= -0.143789 + 0.288874 * v0;
2cf88888 796 Float_t SPDsigma[101]={231.483, 189.446, 183.359, 179.923, 174.229, 170.309, 165.021,
a540a9d3 797 160.84, 159.33, 154.453, 151.644, 148.337, 145.215, 142.353,
798 139.351, 136, 133.838, 129.885, 127.36, 125.032, 122.21, 120.3,
799 117.766, 114.77, 113.1, 110.268, 107.463, 105.293, 102.845,
800 100.835, 98.9632, 97.3287, 93.6887, 92.1066, 89.3224, 87.8382,
801 86.04, 83.6431, 81.9655, 80.0491, 77.8208, 76.4716, 74.2165,
802 72.2752, 70.4875, 68.9414, 66.8622, 65.022, 63.5134, 61.8228,
803 59.7166, 58.5008, 56.2789, 54.9615, 53.386, 51.2165, 49.4842,
804 48.259, 47.1129, 45.3115, 43.8486, 41.9207, 40.5754, 39.3872,
805 38.1897, 36.5401, 35.1283, 33.9702, 32.6429, 31.3612, 29.5876,
806 28.9319, 27.564, 26.0443, 25.2836, 23.9753, 22.8936, 21.5665,
807 20.7048, 19.8016, 18.7095, 18.1144, 17.2095, 16.602, 16.3233,
808 15.7185, 15.3006, 14.7432, 14.4174, 14.0805, 13.7638, 13.7638,
2cf88888 809 13.7638, 13.7638, 13.7638, 13.7638, 13.7638, 13.7638, 13.7638, 18.0803, 18.0803};
a540a9d3 810
f7292781 811 if ( TMath::Abs(spd-val) > fOutliersCut*SPDsigma[cent] )
a540a9d3 812 return kTRUE;
813 else
814 return kFALSE;
815}
816
817//________________________________________________________________________
f690bf48 818Bool_t AliCentralitySelectionTask::IsOutlierV0MTPC(Int_t tracks, Float_t v0, Int_t cent) const
a540a9d3 819{
f690bf48 820// Clean outliers
f7292781 821 Float_t val = -0.540691 + 0.128358 * v0;
2cf88888 822 Float_t TPCsigma[101]={106.439, 89.2834, 86.7568, 85.3641, 83.379, 81.6093, 79.3189,
a540a9d3 823 78.0616, 77.2167, 75.0021, 73.9957, 72.0926, 71.0442, 69.8395,
824 68.1169, 66.6676, 66.0038, 64.2284, 63.3845, 61.7439, 60.642,
825 59.5383, 58.3696, 57.0227, 56.0619, 54.7108, 53.8382, 52.3398,
826 51.5297, 49.9488, 49.1126, 48.208, 46.8566, 45.7724, 44.7829,
827 43.8726, 42.7499, 41.9307, 40.6874, 39.9619, 38.5534, 38.0884,
828 36.6141, 35.7482, 34.8804, 34.1769, 33.1278, 32.3435, 31.4783,
829 30.2587, 29.4741, 28.8575, 27.9298, 26.9752, 26.1675, 25.1234,
830 24.4702, 23.6843, 22.9764, 21.8579, 21.2924, 20.3241, 19.8296,
831 19.2465, 18.4474, 17.7216, 16.8956, 16.342, 15.626, 15.0329,
832 14.3911, 13.9301, 13.254, 12.6745, 12.2436, 11.7776, 11.1795,
833 10.673, 10.27, 9.95646, 9.50939, 9.26162, 8.95315, 8.73439,
834 8.67375, 8.43029, 8.34818, 8.33484, 8.40709, 8.3974, 8.32814,
2cf88888 835 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 12.351, 12.351};
a540a9d3 836
f7292781 837 if ( TMath::Abs(tracks-val) > fOutliersCut*TPCsigma[cent] )
a540a9d3 838 return kTRUE;
839 else
840 return kFALSE;
841}
842
843//________________________________________________________________________
f690bf48 844Bool_t AliCentralitySelectionTask::IsOutlierV0MZDC(Float_t zdc, Float_t v0) const
a540a9d3 845{
f690bf48 846// Clean outliers
f7292781 847 Float_t val1 = 6350. - 0.26 * v0;
848 Float_t val2 = 5580.;
849 if ((zdc > val1) || (zdc > val2))
a540a9d3 850 return kTRUE;
851 else
852 return kFALSE;
853}
854
9138a60c 855//________________________________________________________________________
f690bf48 856Bool_t AliCentralitySelectionTask::IsOutlierV0MZDCECal(Float_t /*zdc*/, Float_t /*v0*/) const
9138a60c 857{
f690bf48 858// Clean outliers
9138a60c 859 return kFALSE;
860}
861
a540a9d3 862//________________________________________________________________________
f690bf48 863Float_t AliCentralitySelectionTask::MyGetScaleFactor(Int_t runnumber, Int_t flag) const
a540a9d3 864{
f690bf48 865// Get scaling factor
a540a9d3 866 if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
867 cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
868 return 0.0;
869 }
870
871 Float_t scalefactor=0.0;
872 if (flag==0)
873 scalefactor = V0MScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
874 else if (flag==1)
875 scalefactor = SPDScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
876 else if (flag==2)
877 scalefactor = TPCScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
878
879 return scalefactor;
880
881}
882
883//________________________________________________________________________
f690bf48 884Float_t AliCentralitySelectionTask::MyGetScaleFactorMC(Int_t runnumber) const
a540a9d3 885{
f690bf48 886// Get MC scaling factor
a540a9d3 887 if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
888 cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
889 return 0.0;
890 }
891
892 Float_t scalefactor= V0MScaleFactorMC[runnumber - fLowRunN]; // subtracting reference offset index
893 return scalefactor;
894
895}
896
897//________________________________________________________________________
898void AliCentralitySelectionTask::MyInitScaleFactor ()
899{
f690bf48 900// Initialize the scaling factors
a540a9d3 901 for (int i=0; i<(fHighRunN-fLowRunN); i++) V0MScaleFactor[i] = 0.0;
902 for (int i=0; i<(fHighRunN-fLowRunN); i++) SPDScaleFactor[i] = 0.0;
903 for (int i=0; i<(fHighRunN-fLowRunN); i++) TPCScaleFactor[i] = 0.0;
904
905 // scale factors determined from <V0 charge> on a run-by-run basis
bbb57131 906 V0MScaleFactor[310] = 1.;
907 V0MScaleFactor[311] = 1.;
609f601e 908 V0MScaleFactor[514] = 1.0046;
909 V0MScaleFactor[515] = 0.983535;
910 V0MScaleFactor[579] = 0.988185;
911 V0MScaleFactor[580] = 0.983351;
912 V0MScaleFactor[581] = 0.989013;
913 V0MScaleFactor[583] = 0.990056;
914 V0MScaleFactor[588] = 0.974438;
915 V0MScaleFactor[589] = 0.981572;
916 V0MScaleFactor[590] = 0.989316;
917 V0MScaleFactor[592] = 0.98345;
918 V0MScaleFactor[688] = 0.993647;
919 V0MScaleFactor[690] = 0.994758;
920 V0MScaleFactor[693] = 0.989569;
921 V0MScaleFactor[698] = 0.993119;
922 V0MScaleFactor[744] = 0.989583;
923 V0MScaleFactor[757] = 0.990377;
924 V0MScaleFactor[787] = 0.990176;
925 V0MScaleFactor[788] = 0.98723;
926 V0MScaleFactor[834] = 1.00403;
927 V0MScaleFactor[835] = 0.994376;
928 V0MScaleFactor[840] = 0.99759;
929 V0MScaleFactor[842] = 1.01031;
930 V0MScaleFactor[900] = 0.996216;
931 V0MScaleFactor[901] = 0.994205;
932 V0MScaleFactor[992] = 0.998479;
933 V0MScaleFactor[997] = 1.00078;
934 V0MScaleFactor[1299] = 1.00515;
935 V0MScaleFactor[1303] = 1.00094;
936 V0MScaleFactor[1339] = 0.986596;
937 V0MScaleFactor[1346] = 0.972226;
938 V0MScaleFactor[1349] = 0.960358;
939 V0MScaleFactor[1350] = 0.970023;
940 V0MScaleFactor[1374] = 1.00575;
941 V0MScaleFactor[1545] = 1.00471;
942 V0MScaleFactor[1587] = 1.00611;
943 V0MScaleFactor[1588] = 1.00976;
944 V0MScaleFactor[1618] = 1.00771;
945 V0MScaleFactor[1682] = 1.01622;
946 V0MScaleFactor[1727] = 1.01305;
947 V0MScaleFactor[1728] = 1.00388;
948 V0MScaleFactor[1731] = 1.00673;
949 V0MScaleFactor[1732] = 1.00916;
950 V0MScaleFactor[1770] = 1.0092;
951 V0MScaleFactor[1773] = 1.00728;
952 V0MScaleFactor[1786] = 1.01655;
953 V0MScaleFactor[1787] = 1.00672;
954 V0MScaleFactor[1801] = 0.983339;
955 V0MScaleFactor[1802] = 1.00754;
956 V0MScaleFactor[1811] = 1.00608;
957 V0MScaleFactor[1815] = 1.01227;
958 V0MScaleFactor[1879] = 0.99907;
959 V0MScaleFactor[1881] = 0.995696;
960 V0MScaleFactor[2186] = 1.00559;
961 V0MScaleFactor[2187] = 1.00631;
962 V0MScaleFactor[2254] = 1.01512;
963 V0MScaleFactor[2256] = 0.998727;
964 V0MScaleFactor[2321] = 1.00701;
965 SPDScaleFactor[310] = 1.00211;
966 SPDScaleFactor[311] = 1.00067;
967 SPDScaleFactor[514] = 1.02268;
968 SPDScaleFactor[515] = 0.994902;
969 SPDScaleFactor[579] = 1.00215;
970 SPDScaleFactor[580] = 0.993421;
971 SPDScaleFactor[581] = 1.00129;
972 SPDScaleFactor[583] = 1.00242;
973 SPDScaleFactor[588] = 0.984762;
974 SPDScaleFactor[589] = 0.994355;
975 SPDScaleFactor[590] = 1.00073;
976 SPDScaleFactor[592] = 0.995889;
977 SPDScaleFactor[688] = 0.994532;
978 SPDScaleFactor[690] = 0.998307;
979 SPDScaleFactor[693] = 0.994052;
980 SPDScaleFactor[698] = 0.993224;
981 SPDScaleFactor[744] = 0.993279;
982 SPDScaleFactor[757] = 0.992494;
983 SPDScaleFactor[787] = 0.992678;
984 SPDScaleFactor[788] = 0.996563;
985 SPDScaleFactor[834] = 1.01116;
986 SPDScaleFactor[835] = 0.993108;
987 SPDScaleFactor[840] = 0.997574;
988 SPDScaleFactor[842] = 1.01829;
989 SPDScaleFactor[900] = 0.999438;
990 SPDScaleFactor[901] = 0.995849;
991 SPDScaleFactor[992] = 0.999227;
992 SPDScaleFactor[997] = 1.00575;
993 SPDScaleFactor[1299] = 0.99877;
994 SPDScaleFactor[1303] = 0.999682;
995 SPDScaleFactor[1339] = 0.978198;
996 SPDScaleFactor[1346] = 0.964178;
997 SPDScaleFactor[1349] = 0.959439;
998 SPDScaleFactor[1350] = 0.956945;
999 SPDScaleFactor[1374] = 0.994434;
1000 SPDScaleFactor[1545] = 1.0016;
1001 SPDScaleFactor[1587] = 1.00153;
1002 SPDScaleFactor[1588] = 1.00698;
1003 SPDScaleFactor[1618] = 1.00554;
1004 SPDScaleFactor[1682] = 1.0123;
1005 SPDScaleFactor[1727] = 1.011;
1006 SPDScaleFactor[1728] = 1.00143;
1007 SPDScaleFactor[1731] = 1.00486;
1008 SPDScaleFactor[1732] = 1.00646;
1009 SPDScaleFactor[1770] = 1.00515;
1010 SPDScaleFactor[1773] = 1.00485;
1011 SPDScaleFactor[1786] = 1.01215;
1012 SPDScaleFactor[1787] = 1.00419;
1013 SPDScaleFactor[1801] = 0.983327;
1014 SPDScaleFactor[1802] = 1.00529;
1015 SPDScaleFactor[1811] = 1.00367;
1016 SPDScaleFactor[1815] = 1.01045;
1017 SPDScaleFactor[1879] = 0.996374;
1018 SPDScaleFactor[1881] = 0.988827;
1019 SPDScaleFactor[2186] = 1.00354;
1020 SPDScaleFactor[2187] = 1.00397;
1021 SPDScaleFactor[2254] = 1.01138;
1022 SPDScaleFactor[2256] = 0.996641;
1023 SPDScaleFactor[2321] = 1.00357;
1024 TPCScaleFactor[310] = 1.00434;
1025 TPCScaleFactor[311] = 1.0056;
1026 TPCScaleFactor[514] = 1.02185;
1027 TPCScaleFactor[515] = 1.0036;
1028 TPCScaleFactor[579] = 1.00607;
1029 TPCScaleFactor[580] = 1.00183;
1030 TPCScaleFactor[581] = 1.00693;
1031 TPCScaleFactor[583] = 1.00746;
1032 TPCScaleFactor[588] = 0.990524;
1033 TPCScaleFactor[589] = 0.998582;
1034 TPCScaleFactor[590] = 1.00618;
1035 TPCScaleFactor[592] = 1.00088;
1036 TPCScaleFactor[688] = 1.00598;
1037 TPCScaleFactor[690] = 1.00658;
1038 TPCScaleFactor[693] = 1.00144;
1039 TPCScaleFactor[698] = 1.00279;
1040 TPCScaleFactor[744] = 1.00122;
1041 TPCScaleFactor[757] = 1.002;
1042 TPCScaleFactor[787] = 0.997818;
1043 TPCScaleFactor[788] = 0.994583;
1044 TPCScaleFactor[834] = 1.01508;
1045 TPCScaleFactor[835] = 1.00218;
1046 TPCScaleFactor[840] = 1.00569;
1047 TPCScaleFactor[842] = 1.01789;
1048 TPCScaleFactor[900] = 1.00739;
1049 TPCScaleFactor[901] = 1.00462;
1050 TPCScaleFactor[992] = 1.00502;
1051 TPCScaleFactor[997] = 1.00943;
1052 TPCScaleFactor[1299] = 1.00438;
1053 TPCScaleFactor[1303] = 0.996701;
1054 TPCScaleFactor[1339] = 0.978641;
1055 TPCScaleFactor[1346] = 0.968906;
1056 TPCScaleFactor[1349] = 0.954311;
1057 TPCScaleFactor[1350] = 0.958764;
1058 TPCScaleFactor[1374] = 0.997899;
1059 TPCScaleFactor[1545] = 0.992;
1060 TPCScaleFactor[1587] = 0.992635;
1061 TPCScaleFactor[1588] = 1.00207;
1062 TPCScaleFactor[1618] = 1.00126;
1063 TPCScaleFactor[1682] = 1.00324;
1064 TPCScaleFactor[1727] = 1.00042;
1065 TPCScaleFactor[1728] = 0.978881;
1066 TPCScaleFactor[1731] = 0.999818;
1067 TPCScaleFactor[1732] = 1.00109;
1068 TPCScaleFactor[1770] = 0.99935;
1069 TPCScaleFactor[1773] = 0.998531;
1070 TPCScaleFactor[1786] = 0.999125;
1071 TPCScaleFactor[1787] = 0.998479;
1072 TPCScaleFactor[1801] = 0.9775;
1073 TPCScaleFactor[1802] = 0.999095;
1074 TPCScaleFactor[1811] = 0.998197;
1075 TPCScaleFactor[1815] = 1.00413;
1076 TPCScaleFactor[1879] = 0.990916;
1077 TPCScaleFactor[1881] = 0.987241;
1078 TPCScaleFactor[2186] = 1.00048;
1079 TPCScaleFactor[2187] = 1.00057;
1080 TPCScaleFactor[2254] = 1.00588;
1081 TPCScaleFactor[2256] = 0.991503;
1082 TPCScaleFactor[2321] = 1.00057;
1083
a540a9d3 1084
1085 // set all missing values to the value of the run before it ....
1086 for (int i=0; i<(fHighRunN-fLowRunN); i++) {
1087 if (V0MScaleFactor[i] == 0.0) {
1088 if (i==0) {
1089 V0MScaleFactor[i] = 1.00;
1090 } else {
1091 // search for last run number with non-zero value
1092 for (int j=i-1;j>=0;j--) {
1093 if (V0MScaleFactor[j] != 0.0) {
1094 V0MScaleFactor[i] = V0MScaleFactor[j];
1095 break;
1096 }
1097 }
1098 }
1099 }
1100 } // end loop over checking all run-numbers
1101
1102 for (int i=0; i<(fHighRunN-fLowRunN); i++) {
1103 if (SPDScaleFactor[i] == 0.0) {
1104 if (i==0) {
1105 SPDScaleFactor[i] = 1.00;
1106 } else {
1107 for (int j=i-1;j>=0;j--) {
1108 if (SPDScaleFactor[j] != 0.0) {
1109 SPDScaleFactor[i] = SPDScaleFactor[j];
1110 break;
1111 }
1112 }
1113 }
1114 }
1115 }
1116
1117 for (int i=0; i<(fHighRunN-fLowRunN); i++) {
1118 if (TPCScaleFactor[i] == 0.0) {
1119 if (i==0) {
1120 TPCScaleFactor[i] = 1.00;
1121 } else {
1122 for (int j=i-1;j>=0;j--) {
1123 if (TPCScaleFactor[j] != 0.0) {
1124 TPCScaleFactor[i] = TPCScaleFactor[j];
1125 break;
1126 }
1127 }
1128 }
1129 }
1130 }
1131
1132
1133 // for (int i=0; i<(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << V0MScaleFactor[i] << " for Run " << i+fLowRunN << endl;
1134
1135 return;
1136
1137}
1138
1139
1140//________________________________________________________________________
1141void AliCentralitySelectionTask::MyInitScaleFactorMC()
1142{
f690bf48 1143// Initialize the MC scaling factors
a540a9d3 1144 for (int i=0; i<(fHighRunN-fLowRunN); i++) V0MScaleFactorMC[i] = 0.0;
1145 // scale factors determined from <V0 charge> on a run-by-run basis
d015c169 1146 V0MScaleFactorMC[0] = 0.75108;
a540a9d3 1147 // set all missing values to the value of the run before it ....
1148 for (int i=0; i<(fHighRunN-fLowRunN); i++) {
1149 if (V0MScaleFactorMC[i] == 0.0) {
1150 if (i==0) {
1151 V0MScaleFactorMC[i] = 1.00;
1152 } else {
1153 // search for last run number with non-zero value
1154 for (int j=i-1;j>=0;j--) {
1155 if (V0MScaleFactorMC[j] != 0.0) {
1156 V0MScaleFactorMC[i] = V0MScaleFactorMC[j];
1157 break;
1158 }
1159 }
1160 }
1161 }
1162 } // end loop over checking all run-numbers
1163
1164
1165 // for (int i=0; i<(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << V0MScaleFactorMC[i] << " for Run " << i+fLowRunN << endl;
1166
1167 return;
1168
1169}
1170