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 //*****************************************************
28 #include <TObjString.h>
32 #include <TDirectory.h>
35 #include "AliAnalysisManager.h"
36 #include "AliVEvent.h"
38 #include "AliESDEvent.h"
39 #include "AliESDHeader.h"
40 #include "AliESDInputHandler.h"
41 #include "AliESDZDC.h"
42 #include "AliESDFMD.h"
43 #include "AliESDVZERO.h"
44 #include "AliESDCentrality.h"
45 #include "AliMultiplicity.h"
46 #include "AliAODHandler.h"
47 #include "AliAODEvent.h"
48 #include "AliAODVertex.h"
49 #include "AliAODMCHeader.h"
50 #include "AliMCEvent.h"
51 #include "AliMCEventHandler.h"
52 #include "AliMCParticle.h"
54 #include "AliHeader.h"
55 #include "AliAODMCParticle.h"
56 #include "AliAnalysisTaskSE.h"
57 #include "AliGenEventHeader.h"
58 #include "AliGenHijingEventHeader.h"
59 #include "AliPhysicsSelectionTask.h"
60 #include "AliPhysicsSelection.h"
61 #include "AliBackgroundSelection.h"
62 #include "AliCentralitySelectionTask.h"
64 ClassImp(AliCentralitySelectionTask)
67 //________________________________________________________________________
68 AliCentralitySelectionTask::AliCentralitySelectionTask():
71 fAnalysisInput("ESD"),
78 fFileList2(new TList),
99 // Default constructor
100 fFileList->SetOwner();
101 fFileList2->SetOwner();
103 AliInfo("Centrality Selection enabled.");
106 //________________________________________________________________________
107 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
108 AliAnalysisTaskSE(name),
110 fAnalysisInput("ESD"),
116 fFileList(new TList),
117 fFileList2(new TList),
138 // Default constructor
139 fFileList->SetOwner();
140 fFileList2->SetOwner();
142 AliInfo("Centrality Selection enabled.");
145 //________________________________________________________________________
146 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
148 // Assignment operator
150 AliAnalysisTaskSE::operator=(c);
155 //________________________________________________________________________
156 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
157 AliAnalysisTaskSE(ana),
159 fAnalysisInput(ana.fDebug),
160 fIsMCInput(ana.fIsMCInput),
163 fCentfilename(ana.fCentfilename),
164 fCentfilename2(ana.fCentfilename2),
165 fFileList(ana.fFileList),
166 fFileList2(ana.fFileList2),
167 fCurrentRun(ana.fCurrentRun),
168 fCentV0M(ana.fCentV0M),
169 fCentFMD(ana.fCentFMD),
170 fCentTRK(ana.fCentTRK),
171 fCentTKL(ana.fCentTKL),
172 fCentCL0(ana.fCentCL0),
173 fCentCL1(ana.fCentCL1),
174 fCentV0MvsFMD(ana.fCentV0MvsFMD),
175 fCentTKLvsV0M(ana.fCentTKLvsV0M),
176 fCentZEMvsZDC(ana.fCentZEMvsZDC),
177 fHtempV0M(ana.fHtempV0M),
178 fHtempFMD(ana.fHtempFMD),
179 fHtempTRK(ana.fHtempTRK),
180 fHtempTKL(ana.fHtempTKL),
181 fHtempCL0(ana.fHtempCL0),
182 fHtempCL1(ana.fHtempCL1),
183 fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
184 fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
185 fHtempZEMvsZDC(ana.fHtempZEMvsZDC)
190 //________________________________________________________________________
191 AliCentralitySelectionTask::~AliCentralitySelectionTask()
208 //________________________________________________________________________
209 void AliCentralitySelectionTask::UserCreateOutputObjects()
211 // Create the output containers
212 if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
214 if (fFileList->GetEntries() < 1) {
215 AliError("No Inputfiles Added");
218 AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
221 //________________________________________________________________________
222 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
224 // Execute analysis for current event:
225 if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
227 Float_t zncEnergy; // ZNC Energy
228 Float_t zpcEnergy; // ZPC Energy
229 Float_t znaEnergy; // ZNA Energy
230 Float_t zpaEnergy; // ZPA Energy
231 Float_t zem1Energy = 0.; // ZEM1 Energy
232 Float_t zem2Energy = 0.; // ZEM2 Energy
234 Int_t nTracks = 0; // no. tracks
235 Int_t nTracklets = 0; // no. tracklets
236 Int_t nClusters[6]; // no. clusters on 6 ITS layers
237 Int_t nChips[2]; // no. chips on 2 SPD layers
239 Float_t multV0A = 0; // multiplicity from V0 reco side A
240 Float_t multV0C = 0; // multiplicity from V0 reco side C
241 Float_t multFMDA = 0; // multiplicity from FMD on detector A
242 Float_t multFMDC = 0; // multiplicity from FMD on detector C
244 AliESDCentrality *esdCent = 0;
246 if(fAnalysisInput.CompareTo("ESD")==0){
248 AliVEvent* event = InputEvent();
249 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
254 esdCent = esd->GetCentrality();
257 AliESDVZERO* esdV0 = esd->GetVZEROData();
258 multV0A=esdV0->GetMTotV0A();
259 multV0C=esdV0->GetMTotV0C();
261 // ***** CB info (tracklets, clusters, chips)
262 nTracks = event->GetNumberOfTracks();
264 const AliMultiplicity *mult = esd->GetMultiplicity();
266 nTracklets = mult->GetNumberOfTracklets();
268 for(Int_t ilay=0; ilay<6; ilay++){
269 nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
272 for(Int_t ilay=0; ilay<2; ilay++){
273 nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
278 AliESDFMD *fmd = esd->GetFMDData();
279 Float_t totalMultA = 0;
280 Float_t totalMultC = 0;
281 const Float_t fFMDLowCut = 0.4;
283 for(UShort_t det=1;det<=3;det++) {
284 Int_t nRings = (det==1 ? 1 : 2);
285 for (UShort_t ir = 0; ir < nRings; ir++) {
286 Char_t ring = (ir == 0 ? 'I' : 'O');
287 UShort_t nsec = (ir == 0 ? 20 : 40);
288 UShort_t nstr = (ir == 0 ? 512 : 256);
289 for(UShort_t sec =0; sec < nsec; sec++) {
290 for(UShort_t strip = 0; strip < nstr; strip++) {
292 Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
293 if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
295 Float_t nParticles=0;
297 if(FMDmult > fFMDLowCut) {
301 if (det<3) totalMultA = totalMultA + nParticles;
302 else totalMultC = totalMultC + nParticles;
308 multFMDA = totalMultA;
309 multFMDC = totalMultC;
312 AliESDZDC *esdZDC = esd->GetESDZDC();
313 zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
314 zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
315 znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
316 zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
317 zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
318 zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
321 else if(fAnalysisInput.CompareTo("AOD")==0){
322 //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
324 printf(" AOD analysis not yet implemented!!!\n\n");
328 // ***** Centrality Selection
329 if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((multV0A+multV0C)));
330 /// else printf(" Centrality by V0 not available!!!\n\n");
331 if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
332 // else printf(" Centrality by FMD not available!!!\n\n");
333 if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
334 // else printf(" Centrality by TRK not available!!!\n\n");
335 if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
336 // else printf(" Centrality by TKL not available!!!\n\n");
337 if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
338 // else printf(" Centrality by CL0 not available!!!\n\n");
339 if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(nClusters[1]));
340 /// else printf(" Centrality by CL1 not available!!!\n\n");
342 if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
343 // else printf(" Centrality by V0 vs FMD not available!!!\n\n");
344 if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
345 // else printf(" Centrality by V0 vs TKL not available!!!\n\n");
346 if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
347 // else printf(" Centrality by ZEM vs ZDC not available!!!\n\n");
349 esdCent->SetCentralityV0M(fCentV0M);
350 esdCent->SetCentralityFMD(fCentFMD);
351 esdCent->SetCentralityTRK(fCentTRK);
352 esdCent->SetCentralityTKL(fCentTKL);
353 esdCent->SetCentralityCL0(fCentCL0);
354 esdCent->SetCentralityCL1(fCentCL1);
355 esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
356 esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
357 esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
360 //________________________________________________________________________
361 void AliCentralitySelectionTask::ReadCentralityHistos()
363 // Read centrality histograms
364 TDirectory *owd = gDirectory;
365 fFile = TFile::Open(fCentfilename);
367 fHtempV0M = (TH1F*) (fFile->Get("hmultV0_percentile"));
368 fHtempFMD = (TH1F*) (fFile->Get("hmultFMD_percentile"));
369 fHtempTRK = (TH1F*) (fFile->Get("hNtracks_percentile"));
370 fHtempTKL = (TH1F*) (fFile->Get("hNtracklets_percentile"));
371 fHtempCL0 = (TH1F*) (fFile->Get("hNclusters0_percentile"));
372 fHtempCL1 = (TH1F*) (fFile->Get("hNclusters1_percentile"));
375 //________________________________________________________________________
376 void AliCentralitySelectionTask::ReadCentralityHistos2()
378 // Read centrality histograms
379 TDirectory *owd = gDirectory;
380 fFile2 = TFile::Open(fCentfilename2);
382 fHtempV0MvsFMD = (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
383 fHtempTKLvsV0M = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
384 fHtempZEMvsZDC = (TH1F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
387 //________________________________________________________________________
388 void AliCentralitySelectionTask::SetPercentileFile(TString filename)
390 // Set the percentile file name
391 fCentfilename = filename;
392 ReadCentralityHistos();
395 //________________________________________________________________________
396 void AliCentralitySelectionTask::SetPercentileFile2(TString filename)
398 // Set the percentile file name
399 fCentfilename2 = filename;
400 ReadCentralityHistos2();
403 //________________________________________________________________________
404 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
406 // Terminate analysis
407 if (fFile && fFile->IsOpen())
409 if (fFile2 && fFile2->IsOpen())
412 //________________________________________________________________________
413 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* esd)
415 // Setup files for run
420 // check if something to be done
421 if (fCurrentRun == esd->GetRunNumber())
424 fCurrentRun = esd->GetRunNumber();
426 AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
428 Int_t runNo = fCurrentRun;
430 // CHANGE HERE FOR RUN RANGES
431 if ( runNo == 137162 ) runNo = 137161;
432 else if ( runNo == 137365 ) runNo = 137366;
433 // CHANGE HERE FOR RUN RANGES
435 TString runName(Form("%d", runNo));
436 TString fileName("");
437 Bool_t isRunKnown = kFALSE;
439 // Check if run is in fileList
440 // if not, take the last name in the list
441 for ( Int_t idx=0 ; idx < fFileList->GetEntries(); ++idx ) {
443 TString str((dynamic_cast<TObjString*>(fFileList->At(idx)))->GetString());
444 if (str.Contains(runName)) {
452 if (fFileList->Last()) {
453 fileName += (dynamic_cast<TObjString*>(fFileList->Last()))->GetString();
454 AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
458 if (fileName.Contains(".root")) {
459 AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));
460 SetPercentileFile(fileName.Data());