1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //*****************************************************
17 // Class AliCentralitySelectionTask
18 // Class to analyze determine centrality
19 // author: Alberica Toia
20 //*****************************************************
22 #include "AliCentralitySelectionTask.h"
31 #include <TObjString.h>
35 #include <TDirectory.h>
39 #include "AliAnalysisManager.h"
40 #include "AliVEvent.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"
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"
70 ClassImp(AliCentralitySelectionTask)
73 //________________________________________________________________________
74 AliCentralitySelectionTask::AliCentralitySelectionTask():
77 fAnalysisInput("ESD"),
115 fHOutCentV0MvsFMD(0),
116 fHOutCentTKLvsV0M(0),
117 fHOutCentZEMvsZDC(0),
125 fHOutMultV0MvsZDC(0),
126 fHOutMultZEMvsZDC(0),
127 fHOutMultV0MvsCL1(0),
128 fHOutMultV0MvsTRK(0),
129 fHOutMultTRKvsCL1(0),
130 fHOutCentV0M_qual1(0),
131 fHOutCentTRK_qual1(0),
132 fHOutCentCL1_qual1(0),
133 fHOutCentV0M_qual2(0),
134 fHOutCentTRK_qual2(0),
135 fHOutCentCL1_qual2(0),
139 // Default constructor
140 AliInfo("Centrality Selection enabled.");
143 for (Int_t i=0; i<(fHighRunN-fLowRunN); i++) {
144 V0MScaleFactor[i]=0.0;
145 SPDScaleFactor[i]=0.0;
146 TPCScaleFactor[i]=0.0;
147 V0MScaleFactorMC[i]=0.0;
152 //________________________________________________________________________
153 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
154 AliAnalysisTaskSE(name),
156 fAnalysisInput("ESD"),
194 fHOutCentV0MvsFMD(0),
195 fHOutCentTKLvsV0M(0),
196 fHOutCentZEMvsZDC(0),
204 fHOutMultV0MvsZDC(0),
205 fHOutMultZEMvsZDC(0),
206 fHOutMultV0MvsCL1(0),
207 fHOutMultV0MvsTRK(0),
208 fHOutMultTRKvsCL1(0),
209 fHOutCentV0M_qual1(0),
210 fHOutCentTRK_qual1(0),
211 fHOutCentCL1_qual1(0),
212 fHOutCentV0M_qual2(0),
213 fHOutCentTRK_qual2(0),
214 fHOutCentCL1_qual2(0),
218 // Default constructor
221 AliInfo("Centrality Selection enabled.");
222 DefineOutput(1, TList::Class());
223 for (Int_t i=0; i<(fHighRunN-fLowRunN); i++) {
224 V0MScaleFactor[i]=0.0;
225 SPDScaleFactor[i]=0.0;
226 TPCScaleFactor[i]=0.0;
227 V0MScaleFactorMC[i]=0.0;
232 //________________________________________________________________________
233 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
235 // Assignment operator
237 AliAnalysisTaskSE::operator=(c);
242 //________________________________________________________________________
243 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
244 AliAnalysisTaskSE(ana),
246 fAnalysisInput(ana.fDebug),
247 fIsMCInput(ana.fIsMCInput),
250 fCurrentRun(ana.fCurrentRun),
252 fLowRunN(ana.fLowRunN),
253 fHighRunN(ana.fHighRunN),
254 fUseScaling(ana.fUseScaling),
255 fTrackCuts(ana.fTrackCuts),
257 fOutliersCut(ana.fOutliersCut),
258 fQuality(ana.fQuality),
259 fCentV0M(ana.fCentV0M),
260 fCentFMD(ana.fCentFMD),
261 fCentTRK(ana.fCentTRK),
262 fCentTKL(ana.fCentTKL),
263 fCentCL0(ana.fCentCL0),
264 fCentCL1(ana.fCentCL1),
265 fCentV0MvsFMD(ana.fCentV0MvsFMD),
266 fCentTKLvsV0M(ana.fCentTKLvsV0M),
267 fCentZEMvsZDC(ana.fCentZEMvsZDC),
268 fHtempV0M(ana.fHtempV0M),
269 fHtempFMD(ana.fHtempFMD),
270 fHtempTRK(ana.fHtempTRK),
271 fHtempTKL(ana.fHtempTKL),
272 fHtempCL0(ana.fHtempCL0),
273 fHtempCL1(ana.fHtempCL1),
274 fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
275 fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
276 fHtempZEMvsZDC(ana.fHtempZEMvsZDC),
277 fOutputList(ana.fOutputList),
278 fHOutCentV0M (ana.fHOutCentV0M ),
279 fHOutCentFMD (ana.fHOutCentFMD ),
280 fHOutCentTRK (ana.fHOutCentTRK ),
281 fHOutCentTKL (ana.fHOutCentTKL ),
282 fHOutCentCL0 (ana.fHOutCentCL0 ),
283 fHOutCentCL1 (ana.fHOutCentCL1 ),
284 fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD),
285 fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M),
286 fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC),
287 fHOutMultV0M(ana.fHOutMultV0M),
288 fHOutMultV0R(ana.fHOutMultV0R),
289 fHOutMultFMD(ana.fHOutMultFMD),
290 fHOutMultTRK(ana.fHOutMultTRK),
291 fHOutMultTKL(ana.fHOutMultTKL),
292 fHOutMultCL0(ana.fHOutMultCL0),
293 fHOutMultCL1(ana.fHOutMultCL1),
294 fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC),
295 fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),
296 fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),
297 fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),
298 fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1),
299 fHOutCentV0M_qual1(ana.fHOutCentV0M_qual1),
300 fHOutCentTRK_qual1(ana.fHOutCentTRK_qual1),
301 fHOutCentCL1_qual1(ana.fHOutCentCL1_qual1),
302 fHOutCentV0M_qual2(ana.fHOutCentV0M_qual2),
303 fHOutCentTRK_qual2(ana.fHOutCentTRK_qual2),
304 fHOutCentCL1_qual2(ana.fHOutCentCL1_qual2),
305 fHOutQuality(ana.fHOutQuality),
306 fHOutVertex(ana.fHOutVertex)
311 //________________________________________________________________________
312 AliCentralitySelectionTask::~AliCentralitySelectionTask()
315 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
316 if (fTrackCuts) delete fTrackCuts;
319 //________________________________________________________________________
320 void AliCentralitySelectionTask::UserCreateOutputObjects()
322 // Create the output containers
323 if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
324 AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
326 fOutputList = new TList();
327 fOutputList->SetOwner();
328 fHOutCentV0M = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",501,0,101);
329 fHOutCentFMD = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",501,0,101);
330 fHOutCentTRK = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",501,0,101);
331 fHOutCentTKL = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",501,0,101);
332 fHOutCentCL0 = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",501,0,101);
333 fHOutCentCL1 = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",501,0,101);
334 fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",501,0,101);
335 fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",501,0,101);
336 fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",501,0,101);
338 fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,25000);
339 fHOutMultV0R = new TH1F("fHOutMultV0R","fHOutMultV0R; Multiplicity V0",25000,0,25000);
340 fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000);
341 fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000);
342 fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000);
343 fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000);
344 fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000);
345 fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,25000,500,0,6000);
346 fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,6000);
347 fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
348 fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
349 fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
351 fHOutCentV0M_qual1 = new TH1F("fHOutCentV0M_qual1","fHOutCentV0M_qual1; Centrality V0",501,0,100);
352 fHOutCentTRK_qual1 = new TH1F("fHOutCentTRK_qual1","fHOutCentTRK_qual1; Centrality TPC",501,0,100);
353 fHOutCentCL1_qual1 = new TH1F("fHOutCentCL1_qual1","fHOutCentCL1_qual1; Centrality SPD outer",501,0,100);
355 fHOutCentV0M_qual2 = new TH1F("fHOutCentV0M_qual2","fHOutCentV0M_qual2; Centrality V0",101,-0.1,100.1);
356 fHOutCentTRK_qual2 = new TH1F("fHOutCentTRK_qual2","fHOutCentTRK_qual2; Centrality TPC",101,-0.1,100.1);
357 fHOutCentCL1_qual2 = new TH1F("fHOutCentCL1_qual2","fHOutCentCL1_qual2; Centrality SPD outer",101,-0.1,100.1);
359 fHOutQuality = new TH1F("fHOutQuality", "fHOutQuality", 10,-0.5,9.5);
360 fHOutVertex = new TH1F("fHOutVertex", "fHOutVertex", 100,-20,20);
362 fOutputList->Add( fHOutCentV0M );
363 fOutputList->Add( fHOutCentFMD );
364 fOutputList->Add( fHOutCentTRK );
365 fOutputList->Add( fHOutCentTKL );
366 fOutputList->Add( fHOutCentCL0 );
367 fOutputList->Add( fHOutCentCL1 );
368 fOutputList->Add( fHOutCentV0MvsFMD);
369 fOutputList->Add( fHOutCentTKLvsV0M);
370 fOutputList->Add( fHOutCentZEMvsZDC);
371 fOutputList->Add( fHOutMultV0M);
372 fOutputList->Add( fHOutMultV0R);
373 fOutputList->Add( fHOutMultFMD);
374 fOutputList->Add( fHOutMultTRK);
375 fOutputList->Add( fHOutMultTKL);
376 fOutputList->Add( fHOutMultCL0);
377 fOutputList->Add( fHOutMultCL1);
378 fOutputList->Add( fHOutMultV0MvsZDC);
379 fOutputList->Add( fHOutMultZEMvsZDC);
380 fOutputList->Add( fHOutMultV0MvsCL1);
381 fOutputList->Add( fHOutMultV0MvsTRK);
382 fOutputList->Add( fHOutMultTRKvsCL1);
383 fOutputList->Add( fHOutCentV0M_qual1 );
384 fOutputList->Add( fHOutCentTRK_qual1 );
385 fOutputList->Add( fHOutCentCL1_qual1 );
386 fOutputList->Add( fHOutCentV0M_qual2 );
387 fOutputList->Add( fHOutCentTRK_qual2 );
388 fOutputList->Add( fHOutCentCL1_qual2 );
389 fOutputList->Add( fHOutQuality );
390 fOutputList->Add( fHOutVertex );
393 fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
395 PostData(1, fOutputList);
398 if (fIsMCInput) MyInitScaleFactorMC();
401 //________________________________________________________________________
402 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
404 // Execute analysis for current event:
405 if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
407 Float_t zncEnergy = 0.; // ZNC Energy
408 Float_t zpcEnergy = 0.; // ZPC Energy
409 Float_t znaEnergy = 0.; // ZNA Energy
410 Float_t zpaEnergy = 0.; // ZPA Energy
411 Float_t zem1Energy = 0.; // ZEM1 Energy
412 Float_t zem2Energy = 0.; // ZEM2 Energy
414 Int_t nTracks = 0; // no. tracks
415 Int_t nTracklets = 0; // no. tracklets
416 Int_t nClusters[6] = {0}; // no. clusters on 6 ITS layers
417 Int_t nChips[2]; // no. chips on 2 SPD layers
418 Float_t spdCorr =0; // corrected spd2 multiplicity
420 Float_t multV0A = 0; // multiplicity from V0 reco side A
421 Float_t multV0C = 0; // multiplicity from V0 reco side C
422 Float_t multFMDA = 0; // multiplicity from FMD on detector A
423 Float_t multFMDC = 0; // multiplicity from FMD on detector C
425 Short_t v0Corr = 0; // corrected V0 multiplicity
426 Short_t v0CorrResc = 0; // corrected and rescaled V0 multiplicity
428 Float_t zvtx =0; // z-vertex SPD
430 AliCentrality *esdCent = 0;
432 if(fAnalysisInput.CompareTo("ESD")==0){
434 AliVEvent* event = InputEvent();
435 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
437 AliError("No ESD Event");
443 AliFatal("Centrality File not available for this run");
446 esdCent = esd->GetCentrality();
449 AliESDVZERO* esdV0 = esd->GetVZEROData();
450 multV0A=esdV0->GetMTotV0A();
451 multV0C=esdV0->GetMTotV0C();
454 v0Corr = (Short_t)AliESDUtils::GetCorrV0(esd,v0CorrR);
455 v0CorrResc = (Short_t)v0CorrR;
458 const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
459 zvtx = vtxESD->GetZ();
461 // ***** CB info (tracklets, clusters, chips)
462 //nTracks = event->GetNumberOfTracks();
463 nTracks = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
465 const AliMultiplicity *mult = esd->GetMultiplicity();
467 nTracklets = mult->GetNumberOfTracklets();
469 for(Int_t ilay=0; ilay<6; ilay++){
470 nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
473 for(Int_t ilay=0; ilay<2; ilay++){
474 nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
477 spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);
480 AliESDFMD *fmd = esd->GetFMDData();
481 Float_t totalMultA = 0;
482 Float_t totalMultC = 0;
483 const Float_t fFMDLowCut = 0.4;
485 for(UShort_t det=1;det<=3;det++) {
486 Int_t nRings = (det==1 ? 1 : 2);
487 for (UShort_t ir = 0; ir < nRings; ir++) {
488 Char_t ring = (ir == 0 ? 'I' : 'O');
489 UShort_t nsec = (ir == 0 ? 20 : 40);
490 UShort_t nstr = (ir == 0 ? 512 : 256);
491 for(UShort_t sec =0; sec < nsec; sec++) {
492 for(UShort_t strip = 0; strip < nstr; strip++) {
494 Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
495 if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
497 Float_t nParticles=0;
499 if(FMDmult > fFMDLowCut) {
503 if (det<3) totalMultA = totalMultA + nParticles;
504 else totalMultC = totalMultC + nParticles;
510 multFMDA = totalMultA;
511 multFMDC = totalMultC;
514 AliESDZDC *esdZDC = esd->GetESDZDC();
515 zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
516 zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
517 znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
518 zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
519 zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
520 zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
523 else if(fAnalysisInput.CompareTo("AOD")==0){
524 //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
526 printf(" AOD analysis not yet implemented!!!\n\n");
533 Float_t temp_scalefactorV0M = MyGetScaleFactor(fCurrentRun,0);
534 Float_t temp_scalefactorSPD = MyGetScaleFactor(fCurrentRun,1);
535 Float_t temp_scalefactorTPC = MyGetScaleFactor(fCurrentRun,2);
536 v0Corr = Short_t(v0Corr / temp_scalefactorV0M);
537 spdCorr = spdCorr / temp_scalefactorSPD;
538 nTracks = Int_t(nTracks / temp_scalefactorTPC);
540 // ***** Scaling for MC
542 Float_t temp_scalefactorV0M = MyGetScaleFactorMC(fCurrentRun);
543 v0Corr = Short_t((multV0A+multV0C) * temp_scalefactorV0M);
546 // ***** Centrality Selection
547 if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
548 if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
549 if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
550 if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
551 if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
552 if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
554 if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
555 if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
556 if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
565 if (TMath::Abs(zvtx)>fZVCut) fQuality += 1;
569 printf("AnalysisCentralitySelectionTask::centrality is %f\n",fCentV0M);
571 if (IsOutlierV0MSPD(spdCorr, v0Corr, int(fCentV0M))) fQuality += 2;
573 if (IsOutlierV0MTPC(nTracks, v0Corr, int(fCentV0M))) fQuality += 4;
575 if (IsOutlierV0MZDC((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr)) fQuality += 8;
578 esdCent->SetQuality(fQuality);
579 esdCent->SetCentralityV0M(fCentV0M);
580 esdCent->SetCentralityFMD(fCentFMD);
581 esdCent->SetCentralityTRK(fCentTRK);
582 esdCent->SetCentralityTKL(fCentTKL);
583 esdCent->SetCentralityCL0(fCentCL0);
584 esdCent->SetCentralityCL1(fCentCL1);
585 esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
586 esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
587 esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
590 fHOutQuality->Fill(fQuality);
591 fHOutVertex->Fill(zvtx);
593 fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
594 fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
595 fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
596 fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
599 fHOutCentV0M->Fill(fCentV0M);
600 fHOutCentFMD->Fill(fCentFMD);
601 fHOutCentTRK->Fill(fCentTRK);
602 fHOutCentTKL->Fill(fCentTKL);
603 fHOutCentCL0->Fill(fCentCL0);
604 fHOutCentCL1->Fill(fCentCL1);
605 fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
606 fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
607 fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
608 fHOutMultV0M->Fill(v0Corr);
609 fHOutMultV0R->Fill(multV0A+multV0C);
610 fHOutMultFMD->Fill((multFMDA+multFMDC));
611 fHOutMultTRK->Fill(nTracks);
612 fHOutMultTKL->Fill(nTracklets);
613 fHOutMultCL0->Fill(nClusters[0]);
614 fHOutMultCL1->Fill(spdCorr);
615 fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
616 } else if (fQuality ==1) {
617 fHOutCentV0M_qual1->Fill(fCentV0M);
618 fHOutCentTRK_qual1->Fill(fCentTRK);
619 fHOutCentCL1_qual1->Fill(fCentCL1);
621 fHOutCentV0M_qual2->Fill(fCentV0M);
622 fHOutCentTRK_qual2->Fill(fCentTRK);
623 fHOutCentCL1_qual2->Fill(fCentCL1);
626 PostData(1, fOutputList);
629 //________________________________________________________________________
630 void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename)
632 // Read centrality histograms
633 TDirectory *owd = gDirectory;
634 // Check if the file is present
635 TString path = gSystem->ExpandPathName(fCentfilename.Data());
636 if (gSystem->AccessPathName(path)) {
637 AliError(Form("File %s does not exist", path.Data()));
640 fFile = TFile::Open(fCentfilename);
642 fHtempV0M = (TH1F*) (fFile->Get("hmultV0_percentile"));
643 fHtempFMD = (TH1F*) (fFile->Get("hmultFMD_percentile"));
644 fHtempTRK = (TH1F*) (fFile->Get("hNtracks_percentile"));
645 fHtempTKL = (TH1F*) (fFile->Get("hNtracklets_percentile"));
646 fHtempCL0 = (TH1F*) (fFile->Get("hNclusters0_percentile"));
647 fHtempCL1 = (TH1F*) (fFile->Get("hNclusters1_percentile"));
649 if (!fHtempV0M) AliWarning(Form("Calibration for V0M does not exist in %s", path.Data()));
650 if (!fHtempFMD) AliWarning(Form("Calibration for FMD does not exist in %s", path.Data()));
651 if (!fHtempTRK) AliWarning(Form("Calibration for TRK does not exist in %s", path.Data()));
652 if (!fHtempTKL) AliWarning(Form("Calibration for TKL does not exist in %s", path.Data()));
653 if (!fHtempCL0) AliWarning(Form("Calibration for CL0 does not exist in %s", path.Data()));
654 if (!fHtempCL1) AliWarning(Form("Calibration for CL1 does not exist in %s", path.Data()));
659 //________________________________________________________________________
660 void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2)
662 // Read centrality histograms
663 TDirectory *owd = gDirectory;
664 TString path = gSystem->ExpandPathName(fCentfilename2.Data());
665 if (gSystem->AccessPathName(path)) {
666 AliError(Form("File %s does not exist", path.Data()));
669 fFile2 = TFile::Open(fCentfilename2);
671 fHtempV0MvsFMD = (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
672 fHtempTKLvsV0M = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
673 fHtempZEMvsZDC = (TH1F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
675 if (!fHtempV0MvsFMD) AliWarning(Form("Calibration for V0MvsFMD does not exist in %s", path.Data()));
676 if (!fHtempTKLvsV0M) AliWarning(Form("Calibration for TKLvsV0M does not exist in %s", path.Data()));
677 if (!fHtempZEMvsZDC) AliWarning(Form("Calibration for ZEMvsZDC does not exist in %s", path.Data()));
682 //________________________________________________________________________
683 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
685 // Terminate analysis
686 if (fFile && fFile->IsOpen())
688 if (fFile2 && fFile2->IsOpen())
691 //________________________________________________________________________
692 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* esd)
694 // Setup files for run
699 // check if something to be done
700 if (fCurrentRun == esd->GetRunNumber())
703 fCurrentRun = esd->GetRunNumber();
705 AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
707 // CHANGE HERE FOR RUN RANGES
708 if ( fCurrentRun <= 137165 ) fRunNo = 137161;
709 else fRunNo = 137366;
710 // CHANGE HERE FOR RUN RANGES
712 TString fileName(Form("%s/COMMON/CENTRALITY/data/AliCentralityBy1D_%d.root", AliAnalysisManager::GetOADBPath(), fRunNo));
713 TString fileName2(Form("%s/COMMON/CENTRALITY/data/AliCentralityByFunction_%d.root", AliAnalysisManager::GetOADBPath(), fRunNo));
715 AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));
716 ReadCentralityHistos(fileName.Data());
717 ReadCentralityHistos2(fileName2.Data());
718 if (!fFile && !fFile2) {
719 AliFatal(Form("Run %d not known to centrality selection!", fCurrentRun));
725 //________________________________________________________________________
726 Bool_t AliCentralitySelectionTask::IsOutlierV0MSPD(Float_t spd, Float_t v0, Int_t cent)
728 TF1 *V0MSPDfun = new TF1("V0MSPDfun","-0.143789+ 0.288874*x",0,25000);
729 Float_t SPDsigma[100]={231.483, 189.446, 183.359, 179.923, 174.229, 170.309, 165.021,
730 160.84, 159.33, 154.453, 151.644, 148.337, 145.215, 142.353,
731 139.351, 136, 133.838, 129.885, 127.36, 125.032, 122.21, 120.3,
732 117.766, 114.77, 113.1, 110.268, 107.463, 105.293, 102.845,
733 100.835, 98.9632, 97.3287, 93.6887, 92.1066, 89.3224, 87.8382,
734 86.04, 83.6431, 81.9655, 80.0491, 77.8208, 76.4716, 74.2165,
735 72.2752, 70.4875, 68.9414, 66.8622, 65.022, 63.5134, 61.8228,
736 59.7166, 58.5008, 56.2789, 54.9615, 53.386, 51.2165, 49.4842,
737 48.259, 47.1129, 45.3115, 43.8486, 41.9207, 40.5754, 39.3872,
738 38.1897, 36.5401, 35.1283, 33.9702, 32.6429, 31.3612, 29.5876,
739 28.9319, 27.564, 26.0443, 25.2836, 23.9753, 22.8936, 21.5665,
740 20.7048, 19.8016, 18.7095, 18.1144, 17.2095, 16.602, 16.3233,
741 15.7185, 15.3006, 14.7432, 14.4174, 14.0805, 13.7638, 13.7638,
742 13.7638, 13.7638, 13.7638, 13.7638, 13.7638, 13.7638, 13.7638, 18.0803};
744 if ( TMath::Abs(spd-V0MSPDfun->Eval(v0)) > fOutliersCut*SPDsigma[cent] )
750 //________________________________________________________________________
751 Bool_t AliCentralitySelectionTask::IsOutlierV0MTPC(Int_t tracks, Float_t v0, Int_t cent)
753 TF1 *V0MTPCfun = new TF1("V0MTPCfun","-0.540691+0.128358*x",0,25000);
754 Float_t TPCsigma[100]={106.439, 89.2834, 86.7568, 85.3641, 83.379, 81.6093, 79.3189,
755 78.0616, 77.2167, 75.0021, 73.9957, 72.0926, 71.0442, 69.8395,
756 68.1169, 66.6676, 66.0038, 64.2284, 63.3845, 61.7439, 60.642,
757 59.5383, 58.3696, 57.0227, 56.0619, 54.7108, 53.8382, 52.3398,
758 51.5297, 49.9488, 49.1126, 48.208, 46.8566, 45.7724, 44.7829,
759 43.8726, 42.7499, 41.9307, 40.6874, 39.9619, 38.5534, 38.0884,
760 36.6141, 35.7482, 34.8804, 34.1769, 33.1278, 32.3435, 31.4783,
761 30.2587, 29.4741, 28.8575, 27.9298, 26.9752, 26.1675, 25.1234,
762 24.4702, 23.6843, 22.9764, 21.8579, 21.2924, 20.3241, 19.8296,
763 19.2465, 18.4474, 17.7216, 16.8956, 16.342, 15.626, 15.0329,
764 14.3911, 13.9301, 13.254, 12.6745, 12.2436, 11.7776, 11.1795,
765 10.673, 10.27, 9.95646, 9.50939, 9.26162, 8.95315, 8.73439,
766 8.67375, 8.43029, 8.34818, 8.33484, 8.40709, 8.3974, 8.32814,
767 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 12.351};
769 if ( TMath::Abs(tracks-V0MTPCfun->Eval(v0)) > fOutliersCut*TPCsigma[cent] )
775 //________________________________________________________________________
776 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDC(Float_t zdc, Float_t v0)
778 TF1 *fun1 = new TF1("fun1","6350-0.26*x",0,25000);
779 TF1 *fun2 = new TF1("fun2","5580",0,25000);
781 if ( (zdc > fun1->Eval(v0)) || (zdc > fun2->Eval(v0)) )
787 //________________________________________________________________________
788 Float_t AliCentralitySelectionTask::MyGetScaleFactor(Int_t runnumber, Int_t flag)
790 if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
791 cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
795 Float_t scalefactor=0.0;
797 scalefactor = V0MScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
799 scalefactor = SPDScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
801 scalefactor = TPCScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
807 //________________________________________________________________________
808 Float_t AliCentralitySelectionTask::MyGetScaleFactorMC(Int_t runnumber)
810 if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
811 cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
815 Float_t scalefactor= V0MScaleFactorMC[runnumber - fLowRunN]; // subtracting reference offset index
820 //________________________________________________________________________
821 void AliCentralitySelectionTask::MyInitScaleFactor ()
823 for (int i=0; i<(fHighRunN-fLowRunN); i++) V0MScaleFactor[i] = 0.0;
824 for (int i=0; i<(fHighRunN-fLowRunN); i++) SPDScaleFactor[i] = 0.0;
825 for (int i=0; i<(fHighRunN-fLowRunN); i++) TPCScaleFactor[i] = 0.0;
827 // scale factors determined from <V0 charge> on a run-by-run basis
828 V0MScaleFactor[0] = 0.956841;
829 V0MScaleFactor[1] = 0.958274;
830 V0MScaleFactor[204] = 1.0046;
831 V0MScaleFactor[205] = 0.983535;
832 V0MScaleFactor[269] = 0.988185;
833 V0MScaleFactor[270] = 0.983351;
834 V0MScaleFactor[271] = 0.989013;
835 V0MScaleFactor[273] = 0.990056;
836 V0MScaleFactor[278] = 0.974438;
837 V0MScaleFactor[279] = 0.981572;
838 V0MScaleFactor[280] = 0.989316;
839 V0MScaleFactor[282] = 0.98345;
840 V0MScaleFactor[378] = 0.993647;
841 V0MScaleFactor[380] = 0.994758;
842 V0MScaleFactor[383] = 0.989569;
843 V0MScaleFactor[388] = 0.993119;
844 V0MScaleFactor[434] = 0.989583;
845 V0MScaleFactor[447] = 0.990377;
846 V0MScaleFactor[477] = 0.990176;
847 V0MScaleFactor[478] = 0.98723;
848 V0MScaleFactor[524] = 1.00403;
849 V0MScaleFactor[525] = 0.994376;
850 V0MScaleFactor[530] = 0.99759;
851 V0MScaleFactor[532] = 1.01031;
852 V0MScaleFactor[590] = 0.996216;
853 V0MScaleFactor[591] = 0.994205;
854 V0MScaleFactor[682] = 0.998479;
855 V0MScaleFactor[687] = 1.00078;
856 V0MScaleFactor[989] = 1.00515;
857 V0MScaleFactor[993] = 1.00094;
858 V0MScaleFactor[1029] = 0.986596;
859 V0MScaleFactor[1036] = 0.972226;
860 V0MScaleFactor[1039] = 0.960358;
861 V0MScaleFactor[1040] = 0.970023;
862 V0MScaleFactor[1064] = 1.00575;
863 V0MScaleFactor[1235] = 1.00471;
864 V0MScaleFactor[1277] = 1.00611;
865 V0MScaleFactor[1278] = 1.00976;
866 V0MScaleFactor[1308] = 1.00771;
867 V0MScaleFactor[1372] = 1.01622;
868 V0MScaleFactor[1417] = 1.01305;
869 V0MScaleFactor[1418] = 1.00388;
870 V0MScaleFactor[1421] = 1.00673;
871 V0MScaleFactor[1422] = 1.00916;
872 V0MScaleFactor[1460] = 1.0092;
873 V0MScaleFactor[1463] = 1.00728;
874 V0MScaleFactor[1476] = 1.01655;
875 V0MScaleFactor[1477] = 1.00672;
876 V0MScaleFactor[1491] = 0.983339;
877 V0MScaleFactor[1492] = 1.00754;
878 V0MScaleFactor[1501] = 1.00608;
879 V0MScaleFactor[1505] = 1.01227;
880 V0MScaleFactor[1569] = 0.99907;
881 V0MScaleFactor[1571] = 0.995696;
882 V0MScaleFactor[1876] = 1.00559;
883 V0MScaleFactor[1877] = 1.00631;
884 V0MScaleFactor[1944] = 1.01512;
885 V0MScaleFactor[1946] = 0.998727;
886 V0MScaleFactor[2011] = 1.00701;
888 SPDScaleFactor[0] = 1.00211;
889 SPDScaleFactor[1] = 1.00067;
890 SPDScaleFactor[204] = 1.02268;
891 SPDScaleFactor[205] = 0.994902;
892 SPDScaleFactor[269] = 1.00215;
893 SPDScaleFactor[270] = 0.993421;
894 SPDScaleFactor[271] = 1.00129;
895 SPDScaleFactor[273] = 1.00242;
896 SPDScaleFactor[278] = 0.984762;
897 SPDScaleFactor[279] = 0.994355;
898 SPDScaleFactor[280] = 1.00073;
899 SPDScaleFactor[282] = 0.995889;
900 SPDScaleFactor[378] = 0.994532;
901 SPDScaleFactor[380] = 0.998307;
902 SPDScaleFactor[383] = 0.994052;
903 SPDScaleFactor[388] = 0.993224;
904 SPDScaleFactor[434] = 0.993279;
905 SPDScaleFactor[447] = 0.992494;
906 SPDScaleFactor[477] = 0.992678;
907 SPDScaleFactor[478] = 0.996563;
908 SPDScaleFactor[524] = 1.01116;
909 SPDScaleFactor[525] = 0.993108;
910 SPDScaleFactor[530] = 0.997574;
911 SPDScaleFactor[532] = 1.01829;
912 SPDScaleFactor[590] = 0.999438;
913 SPDScaleFactor[591] = 0.995849;
914 SPDScaleFactor[682] = 0.999227;
915 SPDScaleFactor[687] = 1.00575;
916 SPDScaleFactor[989] = 0.99877;
917 SPDScaleFactor[993] = 0.999682;
918 SPDScaleFactor[1029] = 0.978198;
919 SPDScaleFactor[1036] = 0.964178;
920 SPDScaleFactor[1039] = 0.959439;
921 SPDScaleFactor[1040] = 0.956945;
922 SPDScaleFactor[1064] = 0.994434;
923 SPDScaleFactor[1235] = 1.0016;
924 SPDScaleFactor[1277] = 1.00153;
925 SPDScaleFactor[1278] = 1.00698;
926 SPDScaleFactor[1308] = 1.00554;
927 SPDScaleFactor[1372] = 1.0123;
928 SPDScaleFactor[1417] = 1.011;
929 SPDScaleFactor[1418] = 1.00143;
930 SPDScaleFactor[1421] = 1.00486;
931 SPDScaleFactor[1422] = 1.00646;
932 SPDScaleFactor[1460] = 1.00515;
933 SPDScaleFactor[1463] = 1.00485;
934 SPDScaleFactor[1476] = 1.01215;
935 SPDScaleFactor[1477] = 1.00419;
936 SPDScaleFactor[1491] = 0.983327;
937 SPDScaleFactor[1492] = 1.00529;
938 SPDScaleFactor[1501] = 1.00367;
939 SPDScaleFactor[1505] = 1.01045;
940 SPDScaleFactor[1569] = 0.996374;
941 SPDScaleFactor[1571] = 0.988827;
942 SPDScaleFactor[1876] = 1.00354;
943 SPDScaleFactor[1877] = 1.00397;
944 SPDScaleFactor[1944] = 1.01138;
945 SPDScaleFactor[1946] = 0.996641;
946 SPDScaleFactor[2011] = 1.00357;
948 TPCScaleFactor[0] = 1.00434;
949 TPCScaleFactor[1] = 1.0056;
950 TPCScaleFactor[204] = 1.02185;
951 TPCScaleFactor[205] = 1.0036;
952 TPCScaleFactor[269] = 1.00607;
953 TPCScaleFactor[270] = 1.00183;
954 TPCScaleFactor[271] = 1.00693;
955 TPCScaleFactor[273] = 1.00746;
956 TPCScaleFactor[278] = 0.990524;
957 TPCScaleFactor[279] = 0.998582;
958 TPCScaleFactor[280] = 1.00618;
959 TPCScaleFactor[282] = 1.00088;
960 TPCScaleFactor[378] = 1.00598;
961 TPCScaleFactor[380] = 1.00658;
962 TPCScaleFactor[383] = 1.00144;
963 TPCScaleFactor[388] = 1.00279;
964 TPCScaleFactor[434] = 1.00122;
965 TPCScaleFactor[447] = 1.002;
966 TPCScaleFactor[477] = 0.997818;
967 TPCScaleFactor[478] = 0.994583;
968 TPCScaleFactor[524] = 1.01508;
969 TPCScaleFactor[525] = 1.00218;
970 TPCScaleFactor[530] = 1.00569;
971 TPCScaleFactor[532] = 1.01789;
972 TPCScaleFactor[590] = 1.00739;
973 TPCScaleFactor[591] = 1.00462;
974 TPCScaleFactor[682] = 1.00502;
975 TPCScaleFactor[687] = 1.00943;
976 TPCScaleFactor[989] = 1.00438;
977 TPCScaleFactor[993] = 0.996701;
978 TPCScaleFactor[1029] = 0.978641;
979 TPCScaleFactor[1036] = 0.968906;
980 TPCScaleFactor[1039] = 0.954311;
981 TPCScaleFactor[1040] = 0.958764;
982 TPCScaleFactor[1064] = 0.997899;
983 TPCScaleFactor[1235] = 0.992;
984 TPCScaleFactor[1277] = 0.992635;
985 TPCScaleFactor[1278] = 1.00207;
986 TPCScaleFactor[1308] = 1.00126;
987 TPCScaleFactor[1372] = 1.00324;
988 TPCScaleFactor[1417] = 1.00042;
989 TPCScaleFactor[1418] = 0.978881;
990 TPCScaleFactor[1421] = 0.999818;
991 TPCScaleFactor[1422] = 1.00109;
992 TPCScaleFactor[1460] = 0.99935;
993 TPCScaleFactor[1463] = 0.998531;
994 TPCScaleFactor[1476] = 0.999125;
995 TPCScaleFactor[1477] = 0.998479;
996 TPCScaleFactor[1491] = 0.9775;
997 TPCScaleFactor[1492] = 0.999095;
998 TPCScaleFactor[1501] = 0.998197;
999 TPCScaleFactor[1505] = 1.00413;
1000 TPCScaleFactor[1569] = 0.990916;
1001 TPCScaleFactor[1571] = 0.987241;
1002 TPCScaleFactor[1876] = 1.00048;
1003 TPCScaleFactor[1877] = 1.00057;
1004 TPCScaleFactor[1944] = 1.00588;
1005 TPCScaleFactor[1946] = 0.991503;
1006 TPCScaleFactor[2011] = 1.00057;
1008 // set all missing values to the value of the run before it ....
1009 for (int i=0; i<(fHighRunN-fLowRunN); i++) {
1010 if (V0MScaleFactor[i] == 0.0) {
1012 V0MScaleFactor[i] = 1.00;
1014 // search for last run number with non-zero value
1015 for (int j=i-1;j>=0;j--) {
1016 if (V0MScaleFactor[j] != 0.0) {
1017 V0MScaleFactor[i] = V0MScaleFactor[j];
1023 } // end loop over checking all run-numbers
1025 for (int i=0; i<(fHighRunN-fLowRunN); i++) {
1026 if (SPDScaleFactor[i] == 0.0) {
1028 SPDScaleFactor[i] = 1.00;
1030 for (int j=i-1;j>=0;j--) {
1031 if (SPDScaleFactor[j] != 0.0) {
1032 SPDScaleFactor[i] = SPDScaleFactor[j];
1040 for (int i=0; i<(fHighRunN-fLowRunN); i++) {
1041 if (TPCScaleFactor[i] == 0.0) {
1043 TPCScaleFactor[i] = 1.00;
1045 for (int j=i-1;j>=0;j--) {
1046 if (TPCScaleFactor[j] != 0.0) {
1047 TPCScaleFactor[i] = TPCScaleFactor[j];
1056 // for (int i=0; i<(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << V0MScaleFactor[i] << " for Run " << i+fLowRunN << endl;
1063 //________________________________________________________________________
1064 void AliCentralitySelectionTask::MyInitScaleFactorMC()
1066 for (int i=0; i<(fHighRunN-fLowRunN); i++) V0MScaleFactorMC[i] = 0.0;
1067 // scale factors determined from <V0 charge> on a run-by-run basis
1068 V0MScaleFactor[0] = 0.75108;
1069 // set all missing values to the value of the run before it ....
1070 for (int i=0; i<(fHighRunN-fLowRunN); i++) {
1071 if (V0MScaleFactorMC[i] == 0.0) {
1073 V0MScaleFactorMC[i] = 1.00;
1075 // search for last run number with non-zero value
1076 for (int j=i-1;j>=0;j--) {
1077 if (V0MScaleFactorMC[j] != 0.0) {
1078 V0MScaleFactorMC[i] = V0MScaleFactorMC[j];
1084 } // end loop over checking all run-numbers
1087 // for (int i=0; i<(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << V0MScaleFactorMC[i] << " for Run " << i+fLowRunN << endl;